Konstancja Satalecka / ctools

Commit 033ed45c7bd6c775c26a8dde33d9d94c973aec44

Authored by Jürgen Knödlseder
2 parents 68b81535 14caa5e9

Merge branch 'integration' into devel

Inline Side-by-side
Showing 60 changed files with 1122 additions and 1814 deletions
pyext/ctool.i
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -25,6 +25,7 @@
25 25 */
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
  28 +#define SWIG
28 29 #include "ctool.hpp"
29 30 %}
30 31  
... ... @@ -51,4 +52,26 @@ public:
51 52  
52 53 // Public methods
53 54 virtual void execute(void);
  55 +
  56 + // Protected methods
  57 + void init_members(void);
  58 + void copy_members(const ctool& app);
  59 + void free_members(void);
  60 + const bool& read_ahead(void) const;
  61 + GObservations get_observations(const bool& get_response = true);
  62 +
  63 + // Protected methods that create objects from user parameters
  64 + GEbounds create_ebounds(void);
  65 + GSkymap create_map(const GObservations& obs);
  66 + GCTAEventCube create_cube(const GObservations& obs);
  67 + GCTAObservation create_cta_obs(void);
  68 +
  69 + // Protected methods that check user parameters
  70 + void require_inobs(const std::string& method);
  71 +
  72 + // Protected support methods
  73 + void set_response(GObservations& obs);
  74 + void set_obs_response(GCTAObservation* obs);
  75 + GSkyDir get_mean_pointing(const GObservations& obs);
  76 + size_t get_current_rss(void);
54 77 };
... ...
pyext/ctools/obsutils.py
  View file @033ed45
... ... @@ -42,6 +42,7 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0,
42 42 binsz - Pixel size for binned simulation (deg/pixel)
43 43 npix - Number of pixels in X and Y for binned simulation
44 44 """
  45 +
45 46 # Allocate ctobssim application and set parameters
46 47 sim = ctools.ctobssim(obs)
47 48 sim["seed"].integer(seed)
... ...
scripts/cspull.py
  View file @033ed45
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # This script generates the pull distribution for all model parameters.
4 4 #
5   -# Copyright (C) 2011-2014 Juergen Knoedlseder
  5 +# Copyright (C) 2011-2015 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -22,7 +22,6 @@ import gammalib
22 22 from ctools import obsutils
23 23 import sys
24 24 import csv
25   -import math
26 25  
27 26  
28 27 # ============ #
... ... @@ -41,12 +40,12 @@ class cspull(gammalib.GApplication):
41 40 """
42 41 # Set name
43 42 self.name = "cspull"
44   - self.version = "0.3.0"
  43 + self.version = "0.4.0"
45 44  
46 45 # Initialise some members
47   - self.obs = None
48   - self.model = None
49   - self.m_srcmdl = None
  46 + self.obs = None
  47 + self.model = None
  48 + self.m_inmodel = None
50 49  
51 50 # Make sure that parfile exists
52 51 file = self.parfile()
... ... @@ -93,7 +92,7 @@ class cspull(gammalib.GApplication):
93 92  
94 93 # Create default parfile
95 94 pars = gammalib.GApplicationPars()
96   - pars.append(gammalib.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
  95 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
97 96 pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
98 97 pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
99 98 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
... ... @@ -123,7 +122,7 @@ class cspull(gammalib.GApplication):
123 122 """
124 123 # Get parameters
125 124 if self.model == None:
126   - self.m_srcmdl = self["srcmdl"].filename()
  125 + self.m_inmodel = self["inmodel"].filename()
127 126 self.m_outfile = self["outfile"].filename()
128 127 self.m_ntrials = self["ntrials"].integer()
129 128 self.m_caldb = self["caldb"].string()
... ... @@ -144,15 +143,15 @@ class cspull(gammalib.GApplication):
144 143  
145 144 # Set some fixed parameters
146 145 self.m_log = False # Logging in client tools
147   - self.m_debug = False # Debugging in client tools
  146 + self.m_debug = True # Debugging in client tools
148 147  
149 148 # Setup observations
150 149 self.obs = self.set_obs()
151   -
  150 +
152 151 # Load source model
153   - if self.m_srcmdl != None:
154   - self.model = gammalib.GModels(self.m_srcmdl)
155   -
  152 + if self.m_inmodel != None:
  153 + self.model = gammalib.GModels(self.m_inmodel)
  154 +
156 155 # Append source model to observation
157 156 self.obs.models(self.model)
158 157  
... ...
scripts/csresmap.py
  View file @033ed45
... ... @@ -21,10 +21,7 @@
21 21 import ctools
22 22 import gammalib
23 23 import sys
24   -import tempfile
25 24 import os
26   -import math
27   -
28 25  
29 26 # ============== #
30 27 # cstsdist class #
... ... @@ -42,7 +39,7 @@ class csresmap(gammalib.GApplication):
42 39 """
43 40 # Set name
44 41 self.name = "csresmap"
45   - self.version = "0.1.0"
  42 + self.version = "0.2.0"
46 43  
47 44 # Initialise some members
48 45 self.obs = None
... ... @@ -103,9 +100,9 @@ class csresmap(gammalib.GApplication):
103 100  
104 101 # Create default parfile
105 102 pars = gammalib.GApplicationPars()
106   - pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition"))
107   - pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name"))
108   - pars.append(gammalib.GApplicationPar("outfile","f","a","resmap.fits","","","Output file name"))
  103 + pars.append(gammalib.GApplicationPar("inobs","f","a","obs.xml","","","Observation definition"))
  104 + pars.append(gammalib.GApplicationPar("inmodel","f","a","crab.xml","","","model file name"))
  105 + pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output file name"))
109 106 pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
110 107 pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
111 108 pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
... ... @@ -137,7 +134,7 @@ class csresmap(gammalib.GApplication):
137 134 # Set observation if not done before
138 135 if self.obs.size() == 0:
139 136  
140   - obsfile = self["infile"].filename()
  137 + obsfile = self["inobs"].filename()
141 138 try:
142 139  
143 140 self.obs = gammalib.GObservations(obsfile)
... ... @@ -163,8 +160,8 @@ class csresmap(gammalib.GApplication):
163 160 # Note that the models should be optimised for
164 161 # meaningful residual map
165 162 if self.obs.models().size() == 0:
166   - self.obs.models(self["srcmdl"].filename())
167   - self.m_outfile = self["outfile"].filename()
  163 + self.obs.models(self["inmodel"].filename())
  164 + self.m_outfile = self["outmap"].filename()
168 165 self.m_modfile = self["modfile"].filename()
169 166 self.m_binfile = self["binfile"].filename()
170 167 self.m_xref = self["xref"].real()
... ...
scripts/cssens.py
  View file @033ed45
... ... @@ -17,7 +17,7 @@
17 17 #
18 18 # The source is modelled as a point source
19 19 #
20   -# Copyright (C) 2011-2014 Juergen Knoedlseder
  20 +# Copyright (C) 2011-2015 Juergen Knoedlseder
21 21 #
22 22 # This program is free software: you can redistribute it and/or modify
23 23 # it under the terms of the GNU General Public License as published by
... ... @@ -33,8 +33,8 @@
33 33 # along with this program. If not, see <http://www.gnu.org/licenses/>.
34 34 #
35 35 # ==========================================================================
36   -from ctools import *
37   -from gammalib import *
  36 +import gammalib
  37 +from ctools import obsutils
38 38 import sys
39 39 import csv
40 40 import math
... ... @@ -43,7 +43,7 @@ import math
43 43 # ============ #
44 44 # cssens class #
45 45 # ============ #
46   -class cssens(GApplication):
  46 +class cssens(gammalib.GApplication):
47 47 """
48 48 This class implements the sensitivity computation script. It derives
49 49 from the GammaLib::GApplication class which provides support for parameter
... ... @@ -63,9 +63,9 @@ class cssens(GApplication):
63 63  
64 64 # Initialise application
65 65 if len(argv) == 0:
66   - GApplication.__init__(self, self.name, self.version)
  66 + gammalib.GApplication.__init__(self, self.name, self.version)
67 67 elif len(argv) ==1:
68   - GApplication.__init__(self, self.name, self.version, *argv)
  68 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
69 69 else:
70 70 raise TypeError("Invalid number of arguments given.")
71 71  
... ... @@ -97,31 +97,31 @@ class cssens(GApplication):
97 97 parfile = self.name+".par"
98 98  
99 99 try:
100   - pars = GApplicationPars(parfile)
  100 + pars = gammalib.GApplicationPars(parfile)
101 101 except:
102 102 # Signal if parfile was not found
103 103 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
104 104  
105 105 # Create default parfile
106   - pars = GApplicationPars()
107   - pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
108   - pars.append(GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database"))
109   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
110   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
111   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
112   - pars.append(GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
113   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
114   - pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
115   - pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
116   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
117   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
118   - pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
119   - pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
120   - pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
121   - pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
122   - pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width"))
123   - pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
124   - pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
  106 + pars = gammalib.GApplicationPars()
  107 + pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
  108 + pars.append(gammalib.GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database"))
  109 + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  110 + pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  111 + pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  112 + pars.append(gammalib.GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  113 + pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  114 + pars.append(gammalib.GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
  115 + pars.append(gammalib.GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
  116 + pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  117 + pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  118 + pars.append(gammalib.GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
  119 + pars.append(gammalib.GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
  120 + pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
  121 + pars.append(gammalib.GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
  122 + pars.append(gammalib.GApplicationPar("width","r","h","0.05","","","Extended source model width"))
  123 + pars.append(gammalib.GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
  124 + pars.append(gammalib.GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
125 125 pars.append_standard()
126 126 pars.save(parfile)
127 127  
... ... @@ -192,8 +192,8 @@ class cssens(GApplication):
192 192 # Initialise models. Note that we centre the point source at the Galactic
193 193 # center as our observation is also centred at the Galactic centre, so
194 194 # we're onaxis.
195   - bkg_model = GModels()
196   - full_model = GModels()
  195 + bkg_model = gammalib.GModels()
  196 + full_model = gammalib.GModels()
197 197 bkg_model.append(self.set_bkg_model())
198 198 full_model.append(self.set_bkg_model())
199 199 full_model.append(self.set_src_model(0.0, \
... ... @@ -263,10 +263,10 @@ class cssens(GApplication):
263 263 emax - Maximum energy [TeV] (default: 100.0)
264 264 """
265 265 # Allocate observation container
266   - obs = GObservations()
  266 + obs = gammalib.GObservations()
267 267  
268 268 # Set single pointing
269   - pntdir = GSkyDir()
  269 + pntdir = gammalib.GSkyDir()
270 270 pntdir.lb_deg(lpnt, bpnt)
271 271  
272 272 # Create CTA observation
... ... @@ -289,21 +289,21 @@ class cssens(GApplication):
289 289 fitidx - Fit spectral index (default: False)
290 290 """
291 291 # Define radial component
292   - radial = GCTAModelRadialGauss(3.0)
  292 + radial = gammalib.GCTAModelRadialGauss(3.0)
293 293  
294 294 # Define spectral component. If a background model is given then
295 295 # use a file function. Otherwise use a static power law.
296 296 if len(self.m_bkg) > 0:
297   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  297 + spectrum = gammalib.GModelSpectralFunc(self.m_bkg, 1.0)
298 298 else:
299   - spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
  299 + spectrum = gammalib.GModelSpectralPlaw(61.8e-6, -1.85, gammalib.GEnergy(1.0, "TeV"))
300 300 if fitidx:
301 301 spectrum["Index"].free()
302 302 else:
303 303 spectrum["Index"].fix()
304 304  
305 305 # Create background model
306   - model = GCTAModelRadialAcceptance(radial, spectrum)
  306 + model = gammalib.GCTAModelRadialAcceptance(radial, spectrum)
307 307 model.name("Background")
308 308 model.instruments("CTA")
309 309  
... ... @@ -328,11 +328,11 @@ class cssens(GApplication):
328 328 fitidx - Fit index? (default: False)
329 329 """
330 330 # Set source location
331   - location = GSkyDir()
  331 + location = gammalib.GSkyDir()
332 332 location.lb_deg(l, b)
333 333  
334 334 # Set source spectrum
335   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
  335 + spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV"))
336 336 if fitidx:
337 337 spectrum["Index"].free()
338 338 else:
... ... @@ -340,7 +340,7 @@ class cssens(GApplication):
340 340  
341 341 # Set source
342 342 if type == "point":
343   - spatial = GModelSpatialPointSource(location)
  343 + spatial = gammalib.GModelSpatialPointSource(location)
344 344 if fitpos:
345 345 spatial["RA"].free()
346 346 spatial["DEC"].free()
... ... @@ -348,7 +348,7 @@ class cssens(GApplication):
348 348 spatial["RA"].fix()
349 349 spatial["DEC"].fix()
350 350 elif type == "gauss":
351   - spatial = GModelSpatialRadialGauss(location, self.m_radius)
  351 + spatial = gammalib.GModelSpatialRadialGauss(location, self.m_radius)
352 352 if fitpos:
353 353 spatial["RA"].free()
354 354 spatial["DEC"].free()
... ... @@ -358,7 +358,7 @@ class cssens(GApplication):
358 358 spatial["DEC"].fix()
359 359 spatial["Sigma"].fix()
360 360 elif type == "disk":
361   - spatial = GModelSpatialRadialDisk(location, self.m_radius)
  361 + spatial = gammalib.GModelSpatialRadialDisk(location, self.m_radius)
362 362 if fitpos:
363 363 spatial["RA"].free()
364 364 spatial["DEC"].free()
... ... @@ -368,7 +368,7 @@ class cssens(GApplication):
368 368 spatial["DEC"].fix()
369 369 spatial["Radius"].fix()
370 370 elif type == "shell":
371   - spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width)
  371 + spatial = gammalib.GModelSpatialRadialShell(location, self.m_radius, self.m_width)
372 372 if fitpos:
373 373 spatial["RA"].free()
374 374 spatial["DEC"].free()
... ... @@ -382,7 +382,7 @@ class cssens(GApplication):
382 382 else:
383 383 self.log("ERROR: Unknown source type '"+type+"'.\n")
384 384 return None
385   - source = GModelSky(spatial, spectrum)
  385 + source = gammalib.GModelSky(spatial, spectrum)
386 386  
387 387 # Set source name
388 388 source.name("Test")
... ...
scripts/cstsdist.py
  View file @033ed45
... ... @@ -3,7 +3,7 @@
3 3 # This script generates the TS distribution for a particular model based
4 4 # on Monte-Carlo simulations.
5 5 #
6   -# Copyright (C) 2011-2014 Juergen Knoedlseder
  6 +# Copyright (C) 2011-2015 Juergen Knoedlseder
7 7 #
8 8 # This program is free software: you can redistribute it and/or modify
9 9 # it under the terms of the GNU General Public License as published by
... ... @@ -19,8 +19,8 @@
19 19 # along with this program. If not, see <http://www.gnu.org/licenses/>.
20 20 #
21 21 # ==========================================================================
22   -from ctools import *
23   -from gammalib import *
  22 +import gammalib
  23 +from ctools import obsutils
24 24 import sys
25 25 import csv
26 26 import math
... ... @@ -29,7 +29,7 @@ import math
29 29 # ============== #
30 30 # cstsdist class #
31 31 # ============== #
32   -class cstsdist(GApplication):
  32 +class cstsdist(gammalib.GApplication):
33 33 """
34 34 This class implements the TS distribution generation script. It derives
35 35 from the GammaLib::GApplication class which provides support for parameter
... ... @@ -42,7 +42,7 @@ class cstsdist(GApplication):
42 42 """
43 43 # Set name
44 44 self.name = "cstsdist"
45   - self.version = "0.1.0"
  45 + self.version = "0.2.0"
46 46  
47 47 # Initialise some members
48 48 self.obs = None
... ... @@ -54,9 +54,9 @@ class cstsdist(GApplication):
54 54  
55 55 # Initialise application
56 56 if len(argv) == 0:
57   - GApplication.__init__(self, self.name, self.version)
  57 + gammalib.GApplication.__init__(self, self.name, self.version)
58 58 elif len(argv) ==1:
59   - GApplication.__init__(self, self.name, self.version, *argv)
  59 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
60 60 else:
61 61 raise TypeError("Invalid number of arguments given.")
62 62  
... ... @@ -87,28 +87,28 @@ class cstsdist(GApplication):
87 87 parfile = self.name+".par"
88 88  
89 89 try:
90   - pars = GApplicationPars(parfile)
  90 + pars = gammalib.GApplicationPars(parfile)
91 91 except:
92 92 # Signal if parfile was not found
93 93 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
94 94  
95 95 # Create default parfile
96   - pars = GApplicationPars()
97   - pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
98   - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
99   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
100   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
101   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
102   - pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index"))
103   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
104   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
105   - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
106   - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
107   - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
108   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
109   - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
110   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
111   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  96 + pars = gammalib.GApplicationPars()
  97 + pars.append(gammalib.GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
  98 + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  99 + pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  100 + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  101 + pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  102 + pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Spectral index"))
  103 + pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  104 + pars.append(gammalib.GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  105 + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  106 + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  107 + pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  108 + pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  109 + pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  110 + pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  111 + pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
112 112 pars.append_standard()
113 113 pars.save(parfile)
114 114  
... ... @@ -146,8 +146,8 @@ class cstsdist(GApplication):
146 146 # Initialise models. Note that we centre the point source at the Galactic
147 147 # center as our observation is also centred at the Galactic centre, so
148 148 # we're onaxis.
149   - self.bkg_model = GModels()
150   - self.full_model = GModels()
  149 + self.bkg_model = gammalib.GModels()
  150 + self.full_model = gammalib.GModels()
151 151 self.bkg_model.append(self.set_bkg_model())
152 152 self.full_model.append(self.set_bkg_model())
153 153 self.full_model.append(self.set_src_model(0.0, self.m_offset, \
... ... @@ -247,10 +247,10 @@ class cstsdist(GApplication):
247 247 emax - Maximum energy [TeV] (default: 100.0)
248 248 """
249 249 # Allocate observation container
250   - obs = GObservations()
  250 + obs = gammalib.GObservations()
251 251  
252 252 # Set single pointing
253   - pntdir = GSkyDir()
  253 + pntdir = gammalib.GSkyDir()
254 254 pntdir.lb_deg(lpnt, bpnt)
255 255  
256 256 # Create CTA observation
... ... @@ -269,17 +269,17 @@ class cstsdist(GApplication):
269 269 Setup CTA background model.
270 270 """
271 271 # Define radial component
272   - radial = GCTAModelRadialGauss(3.0)
  272 + radial = gammalib.GCTAModelRadialGauss(3.0)
273 273 if fitsigma:
274 274 radial["Sigma"].free()
275 275 else:
276 276 radial["Sigma"].fix()
277 277  
278 278 # Define spectral component
279   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  279 + spectrum = gammalib.GModelSpectralFunc(self.m_bkg, 1.0)
280 280  
281 281 # Create background model
282   - model = GCTAModelRadialAcceptance(radial, spectrum)
  282 + model = gammalib.GCTAModelRadialAcceptance(radial, spectrum)
283 283 model.name("Background")
284 284 model.instruments("CTA")
285 285  
... ... @@ -308,11 +308,11 @@ class cstsdist(GApplication):
308 308 fitidx - Fit index? (default: True)
309 309 """
310 310 # Set source location
311   - location = GSkyDir()
  311 + location = gammalib.GSkyDir()
312 312 location.lb_deg(l, b)
313 313  
314 314 # Set source spectrum
315   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
  315 + spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV"))
316 316 if fitidx:
317 317 spectrum["Index"].free()
318 318 else:
... ... @@ -320,32 +320,33 @@ class cstsdist(GApplication):
320 320  
321 321 # Set source
322 322 if type == "point":
323   - spatial = GModelSpatialPointSource(location)
  323 + spatial = gammalib.GModelSpatialPointSource(location)
324 324 if fitpos:
325 325 spatial[0].free()
326 326 spatial[1].free()
327 327 elif type == "gauss":
328   - spatial = GModelSpatialRadialGauss(location, sigma)
  328 + spatial = gammalib.GModelSpatialRadialGauss(location, sigma)
329 329 if fitpos:
330 330 spatial[0].free()
331 331 spatial[1].free()
332 332 elif type == "disk":
333   - spatial = GModelSpatialRadialDisk(location, radius)
  333 + spatial = gammalib.GModelSpatialRadialDisk(location, radius)
334 334 if fitpos:
335 335 spatial[0].free()
336 336 spatial[1].free()
337 337 elif type == "shell":
338   - spatial = GModelSpatialRadialShell(location, radius, width)
  338 + spatial = gammalib.GModelSpatialRadialShell(location, radius, width)
339 339 if fitpos:
340 340 spatial[0].free()
341 341 spatial[1].free()
342 342 else:
343 343 self.log("ERROR: Unknown source type '"+type+"'.\n")
344 344 return None
345   - source = GModelSky(spatial, spectrum)
  345 + source = gammalib.GModelSky(spatial, spectrum)
346 346  
347 347 # Set source name
348 348 source.name("Test")
  349 + source.tscalc(True)
349 350  
350 351 # Return source
351 352 return source
... ... @@ -368,7 +369,7 @@ class cstsdist(GApplication):
368 369 binsz=self.m_binsz, \
369 370 npix=self.m_npix, \
370 371 log=self.m_log, debug=self.m_debug)
371   -
  372 +
372 373 # Determine number of events in simulation
373 374 nevents = 0.0
374 375 for run in sim:
... ... @@ -380,7 +381,7 @@ class cstsdist(GApplication):
380 381 self.log.parformat("Number of simulated events")
381 382 self.log(nevents)
382 383 self.log("\n")
383   -
  384 +
384 385 # Fit background only
385 386 sim.models(self.bkg_model)
386 387 like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
... ...
src/ctbin/ctbin.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctbin - Event binning tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -33,8 +33,8 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36   -#define G_INIT_CUBE "ctbin::init_cube()"
37 36 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)"
  37 +#define G_GET_PARAMETERS "ctbin::get_parameters()"
38 38  
39 39 /* __ Debug definitions __________________________________________________ */
40 40  
... ... @@ -238,9 +238,6 @@ void ctbin::run(void)
238 238 }
239 239 }
240 240  
241   - // Initialise counts cube
242   - init_cube();
243   -
244 241 // Loop over all observations in the container
245 242 for (int i = 0; i < m_obs.size(); ++i) {
246 243  
... ... @@ -309,14 +306,14 @@ void ctbin::save(void)
309 306 }
310 307  
311 308 // Get output filename
312   - m_outfile = (*this)["outfile"].filename();
  309 + m_outcube = (*this)["outcube"].filename();
313 310  
314 311 // Get CTA observation from observation container
315 312 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
316 313  
317 314 // Save only if observation is valid
318 315 if (obs != NULL) {
319   - obs->save(m_outfile, clobber());
  316 + obs->save(m_outcube, clobber());
320 317 }
321 318  
322 319 // Return
... ... @@ -336,16 +333,8 @@ void ctbin::save(void)
336 333 void ctbin::init_members(void)
337 334 {
338 335 // Initialise members
339   - m_evfile.clear();
340   - m_outfile.clear();
341   - m_proj.clear();
342   - m_coordsys.clear();
343   - m_usepnt = false;
344   - m_xref = 0.0;
345   - m_yref = 0.0;
346   - m_binsz = 0.0;
347   - m_nxpix = 0;
348   - m_nypix = 0;
  336 + m_outcube.clear();
  337 + m_usepnt = false;
349 338  
350 339 // Initialise protected members
351 340 m_obs.clear();
... ... @@ -368,24 +357,16 @@ void ctbin::init_members(void)
368 357 void ctbin::copy_members(const ctbin& app)
369 358 {
370 359 // Copy attributes
371   - m_evfile = app.m_evfile;
372   - m_outfile = app.m_outfile;
373   - m_proj = app.m_proj;
374   - m_coordsys = app.m_coordsys;
375   - m_usepnt = app.m_usepnt;
376   - m_xref = app.m_xref;
377   - m_yref = app.m_yref;
378   - m_binsz = app.m_binsz;
379   - m_nxpix = app.m_nxpix;
380   - m_nypix = app.m_nypix;
  360 + m_outcube = app.m_outcube;
  361 + m_usepnt = app.m_usepnt;
381 362  
382 363 // Copy protected members
383   - m_obs = app.m_obs;
384   - m_cube = app.m_cube;
385   - m_ebounds = app.m_ebounds;
386   - m_gti = app.m_gti;
387   - m_ontime = app.m_ontime;
388   - m_livetime = app.m_livetime;
  364 + m_obs = app.m_obs;
  365 + m_cube = app.m_cube;
  366 + m_ebounds = app.m_ebounds;
  367 + m_gti = app.m_gti;
  368 + m_ontime = app.m_ontime;
  369 + m_livetime = app.m_livetime;
389 370  
390 371 // Return
391 372 return;
... ... @@ -413,59 +394,32 @@ void ctbin::free_members(void)
413 394 ***************************************************************************/
414 395 void ctbin::get_parameters(void)
415 396 {
416   - // Initialize energy boundaries
417   - m_ebounds.clear();
418   -
419   - // If there are no observations in container then add a single CTA
420   - // observation using the parameters from the parameter file
  397 + // If there are no observations in container then load them via user
  398 + // parameters
421 399 if (m_obs.size() == 0) {
422 400  
423   - // Get name of CTA events file
424   - m_evfile = (*this)["evfile"].filename();
425   -
426   - // Allocate CTA observation
427   - GCTAObservation obs;
428   -
429   - // Try first to open as FITS file
430   - try {
  401 + // Throw exception if no input observation file is given
  402 + require_inobs(G_GET_PARAMETERS);
431 403  
432   - // Load event list in CTA observation
433   - obs.load(m_evfile);
  404 + // Get observation container without response (not needed)
  405 + m_obs = get_observations(false);
434 406  
435   - // Append CTA observation to container
436   - m_obs.append(obs);
437   -
438   - }
439   -
440   - // ... otherwise try to open as XML file
441   - catch (GException::fits_open_error &e) {
442   -
443   - // Load observations from XML file
444   - m_obs.load(m_evfile);
  407 + } // endif: there was no observation in the container
445 408  
446   - }
  409 + // Create an event cube based on task parameters
  410 + GCTAEventCube cube = create_cube(m_obs);
447 411  
448   - } // endif: there was no observation in the container
  412 + // Get the skymap from the cube and initialise all pixels to zero
  413 + m_cube = cube.map();
  414 + m_cube = 0.0;
449 415  
450   - // Get parameters
451   - m_usepnt = (*this)["usepnt"].boolean();
452   - if (!m_usepnt) {
453   - m_xref = (*this)["xref"].real();
454   - m_yref = (*this)["yref"].real();
455   - }
456   - m_proj = (*this)["proj"].string();
457   - m_coordsys = (*this)["coordsys"].string();
458   - m_binsz = (*this)["binsz"].real();
459   - m_nxpix = (*this)["nxpix"].integer();
460   - m_nypix = (*this)["nypix"].integer();
461   -
462   - // Set energy boundaries
463   - m_ebounds = get_ebounds();
  416 + // Get energy boundaries
  417 + m_ebounds = cube.ebounds();
464 418  
465 419 // Optionally read ahead parameters so that they get correctly
466 420 // dumped into the log file
467 421 if (read_ahead()) {
468   - m_outfile = (*this)["outfile"].filename();
  422 + m_outcube = (*this)["outcube"].filename();
469 423 }
470 424  
471 425 // Return
... ... @@ -474,71 +428,6 @@ void ctbin::get_parameters(void)
474 428  
475 429  
476 430 /***********************************************************************//**
477   - * @brief Initialise counts cube information
478   - *
479   - * @exception GException::invalid_value
480   - * No valid CTA observation found to derive the counts cube
481   - * map centre.
482   - *
483   - * Initialises the skymap, energy boundaries and GTI for a counts cube.
484   - ***************************************************************************/
485   -void ctbin::init_cube(void)
486   -{
487   - // Initialse cube information
488   - m_ontime = 0.0;
489   - m_livetime = 0.0;
490   - m_cube.clear();
491   - m_gti.clear();
492   -
493   - // Set event cube centre, either from the user parameters or from the
494   - // pointing
495   - double xref = m_xref;
496   - double yref = m_yref;
497   - if (m_usepnt) {
498   -
499   - // Dummy: get pointing from first observation. Ultimately, we want
500   - // to get the pointing from a kind of average
501   - bool found = false;
502   - for (int i = 0; i < m_obs.size(); ++i) {
503   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
504   - if (obs != NULL) {
505   - const GCTAPointing& pnt = obs->pointing();
506   - if (gammalib::toupper(m_coordsys) == "GAL") {
507   - xref = pnt.dir().l_deg();
508   - yref = pnt.dir().b_deg();
509   - }
510   - else {
511   - xref = pnt.dir().ra_deg();
512   - yref = pnt.dir().dec_deg();
513   - }
514   - found = true;
515   - break;
516   - }
517   - }
518   -
519   - // Signal if no pointing is found
520   - if (!found) {
521   - std::string msg = "No valid CTA observation has been found in "
522   - "observation list, hence no pointing information "
523   - "could be extracted. Use the \"usepnt=no\" "
524   - "option and specify pointing explicitly when "
525   - "running ctbin.";
526   - throw GException::invalid_value(G_INIT_CUBE, msg);
527   - }
528   -
529   - } // endif: used pointing
530   -
531   - // Create skymap
532   - m_cube = GSkymap(m_proj, m_coordsys,
533   - xref, yref, -m_binsz, m_binsz,
534   - m_nxpix, m_nypix, m_ebounds.size());
535   -
536   - // Return
537   - return;
538   -}
539   -
540   -
541   -/***********************************************************************//**
542 431 * @brief Fill events into counts cube
543 432 *
544 433 * @param[in] obs CTA observation.
... ... @@ -590,8 +479,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
590 479 }
591 480  
592 481 // Skip if pixel is out of range
593   - if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) ||
594   - pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) {
  482 + if (pixel.x() < -0.5 || pixel.x() > (m_cube.nx()-0.5) ||
  483 + pixel.y() < -0.5 || pixel.y() > (m_cube.ny()-0.5)) {
595 484 num_outside_map++;
596 485 continue;
597 486 }
... ...
src/ctbin/ctbin.hpp
  View file @033ed45
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTBIN_NAME "ctbin"
37   -#define CTBIN_VERSION "00-05-00"
  37 +#define CTBIN_VERSION "00-06-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -87,21 +87,12 @@ protected:
87 87 void copy_members(const ctbin& app);
88 88 void free_members(void);
89 89 void get_parameters(void);
90   - void init_cube(void);
91 90 void fill_cube(GCTAObservation* obs);
92 91 void obs_cube(void);
93 92  
94 93 // User parameters
95   - std::string m_evfile; //!< Input event list or XML file
96   - std::string m_outfile; //!< Output counts map or XML file
  94 + std::string m_outcube; //!< Output counts map file
97 95 bool m_usepnt; //!< Use pointing instead of xref/yref parameters
98   - std::string m_proj; //!< WCS projection
99   - std::string m_coordsys; //!< Coordinate system
100   - double m_xref; //!< Longitude reference coordinate
101   - double m_yref; //!< Latitude reference coordinate
102   - double m_binsz; //!< Pixel size
103   - int m_nxpix; //!< Number of pixels in longitude
104   - int m_nypix; //!< Number of pixels in latitude
105 96  
106 97 // Protected members
107 98 GObservations m_obs; //!< Observation container
... ...
src/ctbin/ctbin.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctbin.par #
8   -# Version: 00-05-00 #
  8 +# Version: 00-06-00 #
9 9 # #
10 10 # Author: Jurgen Knodlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -17,8 +17,8 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -evfile, f, a, "events.fits",,,"Input event list or observation definition file"
21   -outfile, f, a, "cntmap.fits",,, "Output counts cube"
  20 +inobs, f, a, "events.fits",,,"Event list or observation definition file"
  21 +outcube, f, a, "cntmap.fits",,, "Output counts cube"
22 22 #
23 23 # Energy binning parameters
24 24 #==========================
... ...
src/ctbkgcube/ctbkgcube.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctbkgcube - Background cube generation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2015 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,8 +34,6 @@
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36 36 #define G_GET_PARAMETERS "ctbkgcube::get_parameters()"
37   -#define G_GET_OBS "ctbkgcube::get_obs()"
38   -#define G_SET_FROM_CNTMAP "ctbkgcube::set_from_cntmap(std::string&)"
39 37 #define G_FILL_CUBE "ctbkgcube::fill_cube(GCTAObservation*)"
40 38  
41 39 /* __ Debug definitions __________________________________________________ */
... ... @@ -258,6 +256,7 @@ void ctbkgcube::run(void)
258 256 // container and put all removed components in the output
259 257 // container
260 258 int num = m_bkgmdl.size();
  259 +
261 260 for (int i = num-1; i >= 0; --i) {
262 261  
263 262 // Flag removal
... ... @@ -407,7 +406,7 @@ void ctbkgcube::save(void)
407 406 }
408 407  
409 408 // Get output filenames
410   - std::string outfile = (*this)["outfile"].filename();
  409 + std::string outfile = (*this)["outcube"].filename();
411 410 std::string outmodel = (*this)["outmodel"].filename();
412 411  
413 412 // Create empty FITS file
... ... @@ -460,6 +459,8 @@ void ctbkgcube::save(void)
460 459 void ctbkgcube::init_members(void)
461 460 {
462 461 // Initialise members
  462 + m_outcube.clear();
  463 + m_outmodel.clear();
463 464 m_obs.clear();
464 465 m_bkgcube.clear();
465 466 m_bkgmdl.clear();
... ... @@ -481,6 +482,8 @@ void ctbkgcube::init_members(void)
481 482 void ctbkgcube::copy_members(const ctbkgcube& app)
482 483 {
483 484 // Copy members
  485 + m_outmodel = app.m_outmodel;
  486 + m_outcube = app.m_outcube;
484 487 m_obs = app.m_obs;
485 488 m_bkgcube = app.m_bkgcube;
486 489 m_bkgmdl = app.m_bkgmdl;
... ... @@ -515,11 +518,18 @@ void ctbkgcube::free_members(void)
515 518 ***************************************************************************/
516 519 void ctbkgcube::get_parameters(void)
517 520 {
518   - // If we do not have any observations in the container then load them
519   - // using the "infile" parameter
  521 + // If there are no observations in container then load them via user
  522 + // parameters
520 523 if (m_obs.size() == 0) {
521   - get_obs();
522   - }
  524 +
  525 + // Throw exception if no input observation file is given
  526 + require_inobs(G_GET_PARAMETERS);
  527 +
  528 + // Build observation container
  529 + m_obs = get_observations();
  530 +
  531 + } // endif: there was no observation in the container
  532 +
523 533  
524 534 // If there are no models associated with the observations then load now
525 535 // the model definition from the XML file
... ... @@ -536,176 +546,47 @@ void ctbkgcube::get_parameters(void)
536 546 m_obs.models(models);
537 547 }
538 548  
539   - // If no counts map is specified then setup the background cube from
540   - // the user parameters
541   - std::string cntmap = (*this)["cntmap"].filename();
542   - if ((gammalib::toupper(cntmap) == "NONE") ||
543   - (gammalib::strip_whitespace(cntmap) == "")) {
544   -
545   - // Get user parameters for counts map definition
546   - std::string proj = (*this)["proj"].string();
547   - std::string coordsys = (*this)["coordsys"].string();
548   - double xref = (*this)["xref"].real();
549   - double yref = (*this)["yref"].real();
550   - double binsz = (*this)["binsz"].real();
551   - int nxpix = (*this)["nxpix"].integer();
552   - int nypix = (*this)["nypix"].integer();
553   -
554   - // Get energy definition
555   - m_ebounds = get_ebounds();
556   -
557   - // Define skymap for background cube
558   - GSkymap map(proj, coordsys, xref, yref,
559   - -binsz, binsz, nxpix, nypix,
560   - m_ebounds.size());
561   -
562   - // Allocate background cube using a dummy GTI
563   - GGti gti;
564   - gti.append(GTime(0.0), GTime(1.0)); // Dummy GTI
565   - m_bkgcube = GCTAEventCube(map, m_ebounds, gti);
  549 + // Get the incube filename
  550 + std::string incube = (*this)["incube"].filename();
566 551  
567   - }
  552 + // Check for filename validity
  553 + if ((gammalib::toupper(incube) == "NONE") ||
  554 + (gammalib::strip_whitespace(incube) == "")) {
  555 +
  556 + // Create an event cube based on task parameters
  557 + m_bkgcube = create_cube(m_obs);
  558 +
  559 + } // endif: filename was not valid
568 560  
569 561 // ... otherwise setup the background cube from the counts map
570 562 else {
571 563  
572   - // Set background cube from counts map
573   - set_from_cntmap(cntmap);
574   -
575   - }
576   -
577   - // Read output filenames (if needed)
578   - if (read_ahead()) {
579   - (*this)["outfile"].filename();
580   - (*this)["outmodel"].filename();
581   - }
582   -
583   - // Return
584   - return;
585   -}
586   -
587   -
588   -/***********************************************************************//**
589   - * @brief Get observation definition
590   - *
591   - * @exception GException::invalid_value
592   - * No event list, counts cube or observation definition file
593   - * specified.
594   - *
595   - * Get observation definition from the user parameters.
596   - ***************************************************************************/
597   -void ctbkgcube::get_obs(void)
598   -{
599   - // Get input filename
600   - std::string filename = (*this)["infile"].filename();
601   -
602   - // Check whether infile parameter is undefined
603   - if (gammalib::toupper(filename) == "NONE") {
604   - std::string msg = "No event list, counts cube or observation definition "
605   - "file specified. Please set the \"infile\" parameter "
606   - "to either an event list filename, a counts cube "
607   - "filename or an observation definition filename.";
608   - throw GException::invalid_value(G_GET_OBS, msg);
609   - }
610   -
611   - // Try first to open as FITS file
612   - try {
613   -
614   - // Allocate CTA observation
615   - GCTAObservation obs;
616   -
617   - // Load input file in CTA observation
618   - obs.load(filename);
619   -
620   - // Append CTA observation to container
621   - m_obs.append(obs);
622   -
623   - }
624   -
625   - // ... otherwise try to open as XML file
626   - catch (GException::fits_open_error &e) {
627   -
628   - // Load observations from XML file
629   - m_obs.load(filename);
630   -
631   - }
632   -
633   - // Return
634   - return;
635   -}
636   -
637   -
638   -/***********************************************************************//**
639   - * @brief Set background cube definition from counts map
640   - *
641   - * @exception GException::invalid_value
642   - * Invalid counts map projection or invalid events encountered.
643   - *
644   - * Set background cube definition from counts map.
645   - ***************************************************************************/
646   -void ctbkgcube::set_from_cntmap(const std::string& filename)
647   -{
648   - // Allocate CTA observation
649   - GCTAObservation obs;
650   -
651   - // Load counts map in CTA observation
652   - obs.load(filename);
653   -
654   - // Set background cube from counts map
655   - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
  564 + // Load event cube from filename
  565 + m_bkgcube.load(incube);
  566 + // Loop over all bins in background cube
  567 + for (int i = 0; i < m_bkgcube.size(); ++i) {
656 568  
657   - // Continue only if cube is valid
658   - if (cube != NULL) {
  569 + // Get event bin
  570 + GCTAEventBin* bin = m_bkgcube[i];
  571 + bin->counts(0.0);
659 572  
660   - // Get sky map projection
661   - const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
662   -
663   - // Continue only if projection is valid
664   - if (wcs != NULL) {
665   -
666   - // Get user parameters for counts map definition
667   - std::string proj = wcs->code();
668   - std::string coordsys = wcs->coordsys();
669   - double xref = wcs->crval(0);
670   - double yref = wcs->crval(1);
671   - double dx = wcs->cdelt(0);
672   - double dy = wcs->cdelt(1);
673   - int nx = cube->map().nx();
674   - int ny = cube->map().ny();
675   -
676   - // Get energy definition
677   - m_ebounds = cube->ebounds();
678   -
679   - // Define skymap for background cube
680   - GSkymap map(proj, coordsys, xref, yref,
681   - dx, dy, nx, ny,
682   - m_ebounds.size());
683   -
684   - // Allocate background cube
685   - m_bkgcube = GCTAEventCube(map, m_ebounds, obs.events()->gti());
686   -
687   - } // endif: WCS projection was valid
688   -
689   - // ... projection is not of WCS type
690   - else {
691   - std::string msg = "Counts map project is not of WCS type.";
692   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
693 573 }
  574 +
  575 + } // endelse: cube was loaded from file
694 576  
695   - } // endif: observation contained an events cube
  577 + // Get energy definition
  578 + m_ebounds = m_bkgcube.ebounds();
696 579  
697   - // ... there is not events cube
698   - else {
699   - std::string msg = "No events cube found in file \""
700   - ""+filename+"\".";
701   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  580 + // Read output filenames (if needed)
  581 + if (read_ahead()) {
  582 + m_outcube = (*this)["outcube"].filename();
  583 + m_outmodel = (*this)["outmodel"].filename();
702 584 }
703 585  
704 586 // Return
705 587 return;
706 588 }
707 589  
708   -
709 590 /***********************************************************************//**
710 591 * @brief Generate background cube
711 592 *
... ... @@ -717,6 +598,7 @@ void ctbkgcube::set_from_cntmap(const std::string&amp; filename)
717 598 ***************************************************************************/
718 599 void ctbkgcube::fill_cube(GCTAObservation* obs)
719 600 {
  601 +
720 602 // Continue only if observation pointer is valid
721 603 if (obs != NULL) {
722 604  
... ...
src/ctbkgcube/ctbkgcube.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctbkgcube - Background cube generation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2015 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTBKGCUBE_NAME "ctbkgcube"
37   -#define CTBKGCUBE_VERSION "00-01-00"
  37 +#define CTBKGCUBE_VERSION "00-02-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -68,11 +68,11 @@ protected:
68 68 void copy_members(const ctbkgcube& app);
69 69 void free_members(void);
70 70 void get_parameters(void);
71   - void get_obs(void);
72   - void set_from_cntmap(const std::string& filename);
73 71 void fill_cube(GCTAObservation* obs);
74 72  
75 73 // Protected members
  74 + std::string m_outcube; //!< Filename of output cube
  75 + std::string m_outmodel; //!< Filename of output XML model
76 76 GObservations m_obs; //!< Observation container
77 77 GCTAEventCube m_bkgcube; //!< Background cube
78 78 GModels m_bkgmdl; //!< CTA background models
... ...
src/ctbkgcube/ctbkgcube.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctbkgcube.par #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Chia-Chun Lu #
11 11 # c.lu@mpi-hd.mpg.de #
... ... @@ -15,13 +15,15 @@
15 15 # #
16 16 #############################################################################
17 17 #
18   -# General parameters
  18 +# Observational parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Event list or observation definition file"
  20 +inobs, f, a, "NONE",,, "Input event list, cube or observation definition file"
21 21 inmodel, f, a, "NONE",,, "Input (background) model XML file"
22   -cntmap, f, a, "NONE",,, "Counts cube for background cube definition"
23   -outfile, f, a, "bkgcube.fits",,, "Output background cube file"
  22 +incube, f, a, "NONE",,, "Counts cube for background cube definition"
  23 +outcube, f, a, "bkgcube.fits",,, "Output background cube file"
24 24 outmodel, f, a, "NONE",,,"Output (background) model XML file"
  25 +caldb, s, a, "dummy",,, "Calibration database"
  26 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 27 #
26 28 # Energy binning parameters
27 29 #==========================
... ... @@ -33,6 +35,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
33 35 #
34 36 # Spatial binning parameters
35 37 #===========================
  38 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
36 39 nxpix, i, a, 200,,, "Size of the X axis in pixels"
37 40 nypix, i, a, 200,,, "Size of the Y axis in pixels"
38 41 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/ctbutterfly/ctbutterfly.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctbutterfly - butterfly calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2014-2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -399,7 +399,6 @@ void ctbutterfly::save(void)
399 399 void ctbutterfly::init_members(void)
400 400 {
401 401 // Initialise members
402   - m_infile.clear();
403 402 m_srcname.clear();
404 403 m_outfile.clear();
405 404 m_ebounds.clear();
... ... @@ -424,7 +423,6 @@ void ctbutterfly::init_members(void)
424 423 void ctbutterfly::copy_members(const ctbutterfly& app)
425 424 {
426 425 // Copy attributes
427   - m_infile = app.m_infile;
428 426 m_srcname = app.m_srcname;
429 427 m_outfile = app.m_outfile;
430 428 m_ebounds = app.m_ebounds;
... ... @@ -466,51 +464,25 @@ void ctbutterfly::free_members(void)
466 464 ***************************************************************************/
467 465 void ctbutterfly::get_parameters(void)
468 466 {
469   - // If there are no observations in container then add a single CTA
470   - // observation using the parameters from the parameter file
  467 + // If there are no observations in container then load them via user
  468 + // parameters
471 469 if (m_obs.size() == 0) {
472 470  
473   - // Allocate CTA observation
474   - GCTAObservation obs;
  471 + // Throw exception if no input observation file is given
  472 + require_inobs(G_GET_PARAMETERS);
475 473  
476   - // Get event file name
477   - std::string filename = (*this)["infile"].filename();
478   -
479   - // Try first to open as FITS file
480   - try {
481   -
482   - // Load data
483   - obs.load(filename);
484   -
485   - // Set response
486   - set_obs_response(&obs);
487   -
488   - // Append observation to container
489   - m_obs.append(obs);
490   -
491   - }
492   -
493   - // ... otherwise try to open as XML file
494   - catch (GException::fits_open_error &e) {
495   -
496   - // Load observations from XML file
497   - m_obs.load(filename);
498   -
499   - // Check if all observations have response information. If
500   - // not, get the calibration database parameters and set
501   - // the response properly
502   - set_response(m_obs);
503   -
504   - } // endcatch: file was an XML file
  474 + // Build observation container
  475 + m_obs = get_observations();
505 476  
506 477 } // endif: there was no observation in the container
507 478  
  479 +
508 480 // If there is are no models associated with the observations then
509 481 // load now the model definition
510 482 if (m_obs.models().size() == 0) {
511 483  
512 484 // Get models XML filename
513   - std::string filename = (*this)["srcmdl"].filename();
  485 + std::string filename = (*this)["inmodel"].filename();
514 486  
515 487 // Setup models for optimizing.
516 488 m_obs.models(GModels(filename));
... ... @@ -526,15 +498,15 @@ void ctbutterfly::get_parameters(void)
526 498 throw GException::invalid_value(G_GET_PARAMETERS, msg);
527 499 }
528 500  
529   - // Get energy binning information
530   - m_ebounds = get_ebounds();
  501 + // Create energy boundaries from user parameters
  502 + m_ebounds = create_ebounds();
531 503  
532 504 // Get matrix file name and load if possible
533 505 std::string matrixfilename = (*this)["matrix"].filename();
534 506 if (matrixfilename != "NONE") {
535 507 std::string msg = "Loading of matrix from file not implemented yet. "
536   - "Use filename = \"NONE\" to induce a recomputation of "
537   - "the matrix internally.";
  508 + "Use filename = \"NONE\" to induce a recomputation of "
  509 + "the matrix internally.";
538 510 throw GException::feature_not_implemented(G_GET_PARAMETERS, msg);
539 511 // m_covariance.load(matrixfilename);
540 512 }
... ...
src/ctbutterfly/ctbutterfly.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctbutterfly - butterfly calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2014-2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTBUTTERFLY_NAME "ctbutterfly"
37   -#define CTBUTTERFLY_VERSION "00-01-00"
  37 +#define CTBUTTERFLY_VERSION "00-02-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -80,12 +80,11 @@ protected:
80 80 void get_parameters(void);
81 81  
82 82 // User parameters
83   - std::string m_infile; //!< Input file
84 83 std::string m_srcname; //!< Name of source to compute butterfly
85 84 std::string m_outfile; //!< Output ascii file
86 85 GEbounds m_ebounds; //!< Energy binning definition
87 86  
88   - // Protected members
  87 + // Protected members
89 88 GObservations m_obs; //!< Observation container
90 89 GMatrixSparse m_covariance; //!< Covariance matrix
91 90 std::vector<double> m_energies; //!< Energy values for storage
... ...
src/ctbutterfly/ctbutterfly.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctbutterfly #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Michael Mayer #
11 11 # michael.mayer@physik.hu-berlin.de #
... ... @@ -17,8 +17,8 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "events.fits",,,"Input event list, cube or observation definition file"
21   -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  20 +inobs, f, a, "events.fits",,,"Input event list, cube or observation definition file"
  21 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
22 22 srcname, s, a, "Crab",,,"Source of interest"
23 23 caldb, s, a, "dummy",,, "Calibration database"
24 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
... ...
src/ctcubemask/ctcubemask.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctcubemask - Cube filter tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2015 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,6 +36,7 @@
36 36 /* __ Method name definitions ____________________________________________ */
37 37 #define G_RUN "ctcubemask::run()"
38 38 #define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)"
  39 +#define G_GET_PARAMETERS "ctcubemask::get_parameters()"
39 40  
40 41 /* __ Debug definitions __________________________________________________ */
41 42  
... ... @@ -336,6 +337,9 @@ void ctcubemask::save(void)
336 337 /***********************************************************************//**
337 338 * @brief Get application parameters
338 339 *
  340 + * @exception GException::invalid_value
  341 + * Parameter "inobs" is required for ctcubemask.
  342 + *
339 343 * Get all task parameters from parameter file or (if required) by querying
340 344 * the user.
341 345 *
... ... @@ -345,46 +349,29 @@ void ctcubemask::save(void)
345 349 ***************************************************************************/
346 350 void ctcubemask::get_parameters(void)
347 351 {
348   - // If there are no observations in container then add a single CTA
349   - // observation using the parameters from the parameter file
  352 + // If there are no observations in container then load them via user
  353 + // parameters
350 354 if (m_obs.size() == 0) {
351 355  
352   - // Get CTA event list file name
353   - m_infile = (*this)["infile"].filename();
354   -
355   - // Allocate CTA observation
356   - GCTAObservation obs;
357   -
358   - // Try first to open as FITS file
359   - try {
360   -
361   - // Load event cube in CTA observation
362   - obs.load(m_infile);
363   -
364   - // Append CTA observation to container
365   - m_obs.append(obs);
366   -
367   - // Signal that no XML file should be used for storage
368   - m_use_xml = false;
369   -
  356 + // Throw exception if no infile is given
  357 + if ((*this)["inobs"].filename() == "NONE" ||
  358 + (*this)["inobs"].filename() == "") {
  359 + std::string msg = "A valid file needs to be specified for the "
  360 + "\"inobs\" parameter, yet \""+
  361 + (*this)["inobs"].filename()+"\" was given."
  362 + " Specify a vaild observation definition or "
  363 + "event list FITS file to proceed.";
  364 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
370 365 }
371   -
372   - // ... otherwise try to open as XML file
373   - catch (GException::fits_open_error &e) {
374   -
375   - // Load observations from XML file
376   - m_obs.load(m_infile);
377 366  
378   - // Signal that XML file should be used for storage
379   - m_use_xml = true;
380   -
381   - }
  367 + // Build observation container without response (not needed)
  368 + m_obs = get_observations(false);
382 369  
383 370 } // endif: there was no observation in the container
384 371  
385 372 // Get parameters
386 373 m_regfile = (*this)["regfile"].filename();
387   - m_usepnt = (*this)["usepnt"].boolean();
  374 + m_usepnt = (*this)["usepnt"].boolean();
388 375 if (!m_usepnt) {
389 376 m_ra = (*this)["ra"].real();
390 377 m_dec = (*this)["dec"].real();
... ... @@ -396,7 +383,7 @@ void ctcubemask::get_parameters(void)
396 383 // Optionally read ahead parameters so that they get correctly
397 384 // dumped into the log file
398 385 if (read_ahead()) {
399   - m_outfile = (*this)["outfile"].filename();
  386 + m_outcube = (*this)["outcube"].filename();
400 387 }
401 388  
402 389 // Return
... ... @@ -538,9 +525,8 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
538 525 void ctcubemask::init_members(void)
539 526 {
540 527 // Initialise parameters
541   - m_infile.clear();
542 528 m_regfile.clear();
543   - m_outfile.clear();
  529 + m_outcube.clear();
544 530 m_prefix.clear();
545 531 m_usepnt = false;
546 532 m_ra = 0.0;
... ... @@ -552,7 +538,6 @@ void ctcubemask::init_members(void)
552 538 // Initialise protected members
553 539 m_obs.clear();
554 540 m_infiles.clear();
555   - m_use_xml = false;
556 541  
557 542 // Return
558 543 return;
... ... @@ -567,9 +552,8 @@ void ctcubemask::init_members(void)
567 552 void ctcubemask::copy_members(const ctcubemask& app)
568 553 {
569 554 // Copy parameters
570   - m_infile = app.m_infile;
571 555 m_regfile = app.m_regfile;
572   - m_outfile = app.m_outfile;
  556 + m_outcube = app.m_outcube;
573 557 m_prefix = app.m_prefix;
574 558 m_usepnt = app.m_usepnt;
575 559 m_ra = app.m_ra;
... ... @@ -581,7 +565,6 @@ void ctcubemask::copy_members(const ctcubemask&amp; app)
581 565 // Copy protected members
582 566 m_obs = app.m_obs;
583 567 m_infiles = app.m_infiles;
584   - m_use_xml = app.m_use_xml;
585 568  
586 569 // Return
587 570 return;
... ... @@ -622,7 +605,7 @@ std::string ctcubemask::check_infile(const std::string&amp; filename) const
622 605 }
623 606 catch (GException::fits_hdu_not_found& e) {
624 607 message = "No \"IMAGE\" extension found in input file \""
625   - + m_outfile + "\".";
  608 + + m_outcube + "\".";
626 609 }
627 610  
628 611 // Return
... ... @@ -661,13 +644,13 @@ std::string ctcubemask::set_outfile_name(const std::string&amp; filename) const
661 644 void ctcubemask::save_fits(void)
662 645 {
663 646 // Get output filename
664   - m_outfile = (*this)["outfile"].filename();
  647 + m_outcube = (*this)["outcube"].filename();
665 648  
666 649 // Get CTA observation from observation container
667 650 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
668 651  
669 652 // Save event list
670   - save_counts_map(obs, m_outfile);
  653 + save_counts_map(obs, m_outcube);
671 654  
672 655 // Return
673 656 return;
... ... @@ -688,14 +671,14 @@ void ctcubemask::save_fits(void)
688 671 void ctcubemask::save_xml(void)
689 672 {
690 673 // Get output filename and prefix
691   - m_outfile = (*this)["outfile"].filename();
  674 + m_outcube = (*this)["m_outcube"].filename();
692 675 m_prefix = (*this)["prefix"].string();
693 676  
694 677 // Issue warning if output filename has no .xml suffix
695   - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
  678 + std::string suffix = gammalib::tolower(m_outcube.substr(m_outcube.length()-4,4));
696 679 if (suffix != ".xml") {
697 680 log << "*** WARNING: Name of observation definition output file \""+
698   - m_outfile+"\"" << std::endl;
  681 + m_outcube+"\"" << std::endl;
699 682 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
700 683 log << "*** WARNING: This is not an error, but might be misleading."
701 684 " It is recommended" << std::endl;
... ... @@ -726,7 +709,7 @@ void ctcubemask::save_xml(void)
726 709 } // endfor: looped over observations
727 710  
728 711 // Save observations in XML file
729   - m_obs.save(m_outfile);
  712 + m_obs.save(m_outcube);
730 713  
731 714 // Return
732 715 return;
... ...
src/ctcubemask/ctcubemask.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctcubemask - Cube filter tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2015 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,7 +36,7 @@
36 36  
37 37 /* __Definitions _________________________________________________________ */
38 38 #define CTCUBEMASK_NAME "ctcubemask"
39   -#define CTCUBEMASK_VERSION "00-01-00"
  39 +#define CTCUBEMASK_VERSION "00-02-00"
40 40  
41 41  
42 42 /***********************************************************************//**
... ... @@ -78,9 +78,9 @@ protected:
78 78 const std::string& outfile) const;
79 79  
80 80 // User parameters
81   - std::string m_infile; //!< Input event list or XML file
  81 + //std::string m_infile; //!< Input event list or XML file
82 82 std::string m_regfile; //!< ds9 region file
83   - std::string m_outfile; //!< Output event list or XML file
  83 + std::string m_outcube; //!< Output event list or XML file
84 84 std::string m_prefix; //!< Prefix for multiple counts maps
85 85 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
86 86 double m_ra; //!< RA of ROI centre
... ... @@ -92,7 +92,6 @@ protected:
92 92 // Protected members
93 93 GObservations m_obs; //!< Observations container
94 94 std::vector<std::string> m_infiles; //!< Input event filenames
95   - bool m_use_xml; //!< Use XML file instead of FITS file
96 95 };
97 96  
98 97  
... ...
src/ctcubemask/ctcubemask.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctcubemask.par #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Chia-Chun Lu #
11 11 # c.lu@mpi-hd.mpg.de #
... ... @@ -17,9 +17,9 @@
17 17 #
18 18 # File information
19 19 #=================
20   -infile, f, a, "cntmap.fits",,,"Input counts cube or observation definition file"
  20 +inobs, f, a, "cntmap.fits",,,"Input counts cube or observation definition file"
21 21 regfile, f, a, "NONE",,,"Exclusion region file in ds9 format"
22   -outfile, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file"
  22 +outcube, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file"
23 23 prefix, s, h, "filtered_",,,"Prefix for counts cube in observation definition file"
24 24 #
25 25 # Selection parameters
... ...
src/ctexpcube/ctexpcube.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctexpcube - Exposure cube generation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -284,10 +284,10 @@ void ctexpcube::save(void)
284 284 }
285 285  
286 286 // Get output filename
287   - m_outfile = (*this)["outfile"].filename();
  287 + m_outcube = (*this)["outcube"].filename();
288 288  
289 289 // Save exposure cube
290   - m_expcube.save(m_outfile, clobber());
  290 + m_expcube.save(m_outcube, clobber());
291 291  
292 292 // Return
293 293 return;
... ... @@ -306,13 +306,12 @@ void ctexpcube::save(void)
306 306 void ctexpcube::init_members(void)
307 307 {
308 308 // Initialise members
309   - m_outfile.clear();
  309 + m_outcube.clear();
310 310 m_apply_edisp = false;
311 311  
312 312 // Initialise protected members
313 313 m_obs.clear();
314 314 m_expcube.clear();
315   - m_ebounds.clear();
316 315  
317 316 // Return
318 317 return;
... ... @@ -327,13 +326,12 @@ void ctexpcube::init_members(void)
327 326 void ctexpcube::copy_members(const ctexpcube& app)
328 327 {
329 328 // Copy attributes
330   - m_outfile = app.m_outfile;
  329 + m_outcube = app.m_outcube;
331 330 m_apply_edisp = app.m_apply_edisp;
332 331  
333 332 // Copy protected members
334   - m_obs = app.m_obs;
335   - m_expcube = app.m_expcube;
336   - m_ebounds = app.m_ebounds;
  333 + m_obs = app.m_obs;
  334 + m_expcube = app.m_expcube;
337 335  
338 336 // Return
339 337 return;
... ... @@ -355,164 +353,52 @@ void ctexpcube::free_members(void)
355 353 *
356 354 * Get all task parameters from parameter file or (if required) by querying
357 355 * the user. The parameters are read in the correct order.
358   - *
359   - * @todo Setup exposure cube from counts map
360 356 ***************************************************************************/
361 357 void ctexpcube::get_parameters(void)
362 358 {
363   - // If we do not have any observations in the container then get an
364   - // input file name or observation descriptor file
  359 + // If there are no observations in container then load them via user
  360 + // parameters
365 361 if (m_obs.size() == 0) {
366   - get_obs();
367   - }
368   -
369   - // For all observations that have no response, set the response
370   - // from the task parameters
371   - set_response(m_obs);
372   -
373   - // Read energy dispersion flag
374   - m_apply_edisp = (*this)["edisp"].boolean();
375 362  
376   - // If no counts map is specified then setup the exposure cube from
377   - // the user parameters
378   - std::string cntmap = (*this)["cntmap"].filename();
379   - if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) {
380   -
381   - // Get user parameters for counts map definition
382   - std::string wcs = (*this)["proj"].string();
383   - std::string coordsys = (*this)["coordsys"].string();
384   - double xref = (*this)["xref"].real();
385   - double yref = (*this)["yref"].real();
386   - double binsz = (*this)["binsz"].real();
387   - int nxpix = (*this)["nxpix"].integer();
388   - int nypix = (*this)["nypix"].integer();
389   -
390   - // Get energy definition
391   - m_ebounds = get_ebounds();
  363 + // Throw exception if no input observation file is given
  364 + require_inobs(G_GET_PARAMETERS);
392 365  
393   - // Define exposure cube
394   - m_expcube = GCTAExposure(wcs, coordsys, xref, yref,
395   - -binsz, binsz, nxpix, nypix,
396   - m_ebounds);
397   - }
  366 + // Build observation container
  367 + m_obs = get_observations();
398 368  
399   - // ... otherwise setup the exposure cube from the counts map
400   - else {
401   -
402   - // Set exposure cube from counts map
403   - set_from_cntmap(cntmap);
404   -
405   - }
  369 + } // endif: there was no observation in the container
406 370  
407   - // Read output filename (if needed)
408   - if (read_ahead()) {
409   - m_outfile = (*this)["outfile"].filename();
410   - }
411   -
412   - // Return
413   - return;
414   -}
415   -
416   -
417   -/***********************************************************************//**
418   - * @brief Get observation definition
419   - *
420   - * Get observation definition from the user parameters.
421   - ***************************************************************************/
422   -void ctexpcube::get_obs(void)
423   -{
424   - // Get input filename
425   - std::string filename = (*this)["infile"].filename();
426   -
427   - // Try first to open as FITS file
428   - try {
  371 + // Read energy dispersion flag
  372 + m_apply_edisp = (*this)["edisp"].boolean();
429 373  
430   - // Allocate CTA observation
431   - GCTAObservation obs;
432   -
433   - // Load input file in CTA observation
434   - obs.load(filename);
  374 + // Get the incube filename
  375 + std::string incube = (*this)["incube"].filename();
435 376  
436   - // Append CTA observation to container
437   - m_obs.append(obs);
  377 + // Check for filename validity
  378 + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
438 379  
439   -
440   - }
441   -
442   - // ... otherwise try to open as XML file
443   - catch (GException::fits_open_error &e) {
  380 + // Create an event cube based on task parameters
  381 + GCTAEventCube cube = create_cube(m_obs);
444 382  
445   - // Load observations from XML file
446   - m_obs.load(filename);
  383 + // Define exposure cube
  384 + m_expcube = GCTAExposure(cube);
447 385  
448   - }
  386 + } // endif: filename was not valid
449 387  
450   - // Return
451   - return;
452   -}
  388 + // ... otherwise setup the exposure cube from the counts map
  389 + else {
453 390  
  391 + // Load event cube from filename
  392 + GCTAEventCube cube(incube);
454 393  
455   -/***********************************************************************//**
456   - * @brief Set exposure cube definition from counts map
457   - *
458   - * @exception GException::invalid_value
459   - * Invalid counts map projection or invalid events encountered.
460   - *
461   - * Set exposure cube definition from counts map.
462   - ***************************************************************************/
463   -void ctexpcube::set_from_cntmap(const std::string& filename)
464   -{
465   - // Allocate CTA observation
466   - GCTAObservation obs;
467   -
468   - // Load counts map in CTA observation
469   - obs.load(filename);
470   -
471   - // Set exposure cube from counts map
472   - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
473   -
474   - // Continue only if cube is valid
475   - if (cube != NULL) {
476   -
477   - // Get sky map projection
478   - const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
479   -
480   - // Continue only if projection is valid
481   - if (wcs != NULL) {
482   -
483   - // Get user parameters for counts map definition
484   - std::string proj = wcs->code();
485   - std::string coordsys = wcs->coordsys();
486   - double xref = wcs->crval(0);
487   - double yref = wcs->crval(1);
488   - double dx = wcs->cdelt(0);
489   - double dy = wcs->cdelt(1);
490   - int nx = cube->map().nx();
491   - int ny = cube->map().ny();
492   -
493   - // Get energy definition
494   - m_ebounds = cube->ebounds();
495   -
496   - // Define exposure cube
497   - m_expcube = GCTAExposure(proj, coordsys, xref, yref,
498   - dx, dy, nx, ny,
499   - m_ebounds);
500   -
501   - } // endif: WCS projection was valid
502   -
503   - // ... projection is not of WCS type
504   - else {
505   - std::string msg = "Counts map project is not of WCS type.";
506   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
507   - }
  394 + // Define exposure cube
  395 + m_expcube = GCTAExposure(cube);
508 396  
509   - } // endif: observation contained an events cube
  397 + } // endelse: cube loaded from file
510 398  
511   - // ... there is not events cube
512   - else {
513   - std::string msg = "No events cube found in file \""
514   - ""+filename+"\".";
515   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  399 + // Read output filename (if needed)
  400 + if (read_ahead()) {
  401 + m_outcube = (*this)["outcube"].filename();
516 402 }
517 403  
518 404 // Return
... ...
src/ctexpcube/ctexpcube.hpp
  View file @033ed45
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTEXPCUBE_NAME "ctexpcube"
37   -#define CTEXPCUBE_VERSION "00-01-00"
  37 +#define CTEXPCUBE_VERSION "00-02-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -67,17 +67,14 @@ protected:
67 67 void copy_members(const ctexpcube& app);
68 68 void free_members(void);
69 69 void get_parameters(void);
70   - void get_obs(void);
71   - void set_from_cntmap(const std::string& filename);
72 70  
73 71 // User parameters
74   - std::string m_outfile; //!< Output exposure cube file
  72 + std::string m_outcube; //!< Output exposure cube file
75 73 bool m_apply_edisp; //!< Apply energy dispersion?
76 74  
77 75 // Protected members
78 76 GObservations m_obs; //!< Observation container
79 77 GCTAExposure m_expcube; //!< Exposure cube
80   - GEbounds m_ebounds; //!< Energy boundaries
81 78 };
82 79  
83 80  
... ...
src/ctexpcube/ctexpcube.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctexpcube.par #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -15,11 +15,11 @@
15 15 # #
16 16 #############################################################################
17 17 #
18   -# General parameters
  18 +# Observational parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Event list or observation definition file"
21   -cntmap, f, a, "NONE",,, "Counts map for exposure cube definition"
22   -outfile, f, a, "expcube.fits",,, "Output exposure cube file"
  20 +inobs, f, a, "NONE",,, "Event list or observation definition file"
  21 +incube, f, a, "NONE",,, "Counts map for exposure cube definition"
  22 +outcube, f, a, "expcube.fits",,, "Output exposure cube file"
23 23 caldb, s, a, "dummy",,, "Calibration database"
24 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 25 edisp, b, h, no,,, "Apply energy dispersion?"
... ... @@ -34,6 +34,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
34 34 #
35 35 # Spatial binning parameters
36 36 #===========================
  37 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
37 38 nxpix, i, a, 200,,, "Size of the X axis in pixels"
38 39 nypix, i, a, 200,,, "Size of the Y axis in pixels"
39 40 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/ctlike/ctlike.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctlike - Maximum likelihood fitting tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -326,11 +326,11 @@ void ctlike::save(void)
326 326 }
327 327  
328 328 // Get output filename
329   - m_outmdl = (*this)["outmdl"].filename();
  329 + m_outmodel = (*this)["outmodel"].filename();
330 330  
331 331 // Write results out as XML model
332   - if (gammalib::toupper(m_outmdl) != "NONE") {
333   - m_obs.models().save(m_outmdl);
  332 + if (gammalib::toupper(m_outmodel) != "NONE") {
  333 + m_obs.models().save(m_outmodel);
334 334 }
335 335  
336 336 // Return
... ... @@ -357,56 +357,34 @@ void ctlike::save(void)
357 357 ***************************************************************************/
358 358 void ctlike::get_parameters(void)
359 359 {
360   - // If there are no observations in container then add a single CTA
361   - // observation using the parameters from the parameter file
  360 + // If there are no observations in container then load them via user
  361 + // parameters
362 362 if (m_obs.size() == 0) {
363 363  
364   - // Allocate CTA observation
365   - GCTAObservation obs;
  364 + // Throw exception if no input observation file is given
  365 + require_inobs(G_GET_PARAMETERS);
366 366  
367   - // Get event file name
368   - std::string filename = (*this)["infile"].filename();
  367 + // Get observation container
  368 + m_obs = get_observations();
369 369  
370   - // Try first to open as FITS file
371   - try {
372   -
373   - // Load data
374   - obs.load(filename);
375   -
376   - // Get other task parameters
377   - std::string statistics = gammalib::toupper((*this)["stat"].string());
378   -
379   - // Set statistics
380   - obs.statistics(statistics);
381   -
382   - // Set response
383   - set_obs_response(&obs);
384   -
385   - // Append observation to container
386   - m_obs.append(obs);
387   -
388   - }
389   -
390   - // ... otherwise try to open as XML file
391   - catch (GException::fits_open_error &e) {
392   -
393   - // Load observations from XML file
394   - m_obs.load(filename);
  370 + } // endif: there was no observation in the container
395 371  
396   - // For all observations that have no response, set the response
397   - // from the task parameters
398   - set_response(m_obs);
  372 + // If only single observation is used, read statistics parameter
  373 + if (!m_use_xml) {
399 374  
400   - } // endcatch: file was an XML file
  375 + // Get other task parameters
  376 + std::string statistics = gammalib::toupper((*this)["stat"].string());
401 377  
402   - } // endif: there was no observation in the container
  378 + // Set statistics
  379 + (*m_obs[0]).statistics(statistics);
  380 + }
403 381  
404 382 // If there is are no models associated with the observations then
405 383 // load now the model definition
406 384 if (m_obs.models().size() == 0) {
407 385  
408 386 // Get models XML filename
409   - std::string filename = (*this)["srcmdl"].filename();
  387 + std::string filename = (*this)["inmodel"].filename();
410 388  
411 389 // Setup models for optimizing.
412 390 m_obs.models(GModels(filename));
... ... @@ -420,7 +398,7 @@ void ctlike::get_parameters(void)
420 398 // Optionally read ahead parameters so that they get correctly
421 399 // dumped into the log file
422 400 if (read_ahead()) {
423   - m_outmdl = (*this)["outmdl"].filename();
  401 + m_outmodel = (*this)["outmodel"].filename();
424 402 }
425 403  
426 404 // Return
... ... @@ -547,7 +525,7 @@ double ctlike::reoptimize_lm(void)
547 525 void ctlike::init_members(void)
548 526 {
549 527 // Initialise members
550   - m_outmdl.clear();
  528 + m_outmodel.clear();
551 529 m_obs.clear();
552 530 m_refit = false;
553 531 m_max_iter = 100; // Set maximum number of iterations
... ... @@ -573,7 +551,7 @@ void ctlike::copy_members(const ctlike&amp; app)
573 551 {
574 552 // Copy attributes
575 553 m_refit = app.m_refit;
576   - m_outmdl = app.m_outmdl;
  554 + m_outmodel = app.m_outmodel;
577 555 m_obs = app.m_obs;
578 556 m_max_iter = app.m_max_iter;
579 557 m_max_stall = app.m_max_stall;
... ...
src/ctlike/ctlike.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctlike - Maximum likelihood fitting tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,7 +35,7 @@
35 35  
36 36 /* __Definitions _________________________________________________________ */
37 37 #define CTLIKE_NAME "ctlike"
38   -#define CTLIKE_VERSION "00-06-00"
  38 +#define CTLIKE_VERSION "00-07-00"
39 39  
40 40  
41 41 /***********************************************************************//**
... ... @@ -74,7 +74,7 @@ protected:
74 74  
75 75 // User parameters
76 76 bool m_refit; //!< Refitting
77   - std::string m_outmdl; //!< Source model output XML file
  77 + std::string m_outmodel; //!< Source model output XML file
78 78 bool m_apply_edisp; //!< Apply energy dispersion?
79 79  
80 80 // Members
... ...
src/ctlike/ctlike.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctlike.par #
8   -# Version: 00-05-00 #
  8 +# Version: 00-06-00 #
9 9 # #
10 10 # Author: Jurgen Knodlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -17,14 +17,14 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "events.fits",,, "Event list, counts cube or observation definition file"
21   -stat, s, h, "POISSON",,, "Optimization statistics"
22   -edisp, b, h, no,,, "Apply energy dispersion?"
23   -refit, b, h, no,,, "Do refitting?"
24   -caldb, s, a, "dummy",,, "Calibration database"
25   -irf, s, a, "cta_dummy_irf",,, "Instrument response function"
26   -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
27   -outmdl, f, a, "crab_results.xml",,, "Source model output file"
  20 +inobs, f, a, "events.fits",,, "Event list, counts cube or observation definition file"
  21 +stat, s, h, "POISSON",,, "Optimization statistics"
  22 +edisp, b, h, no,,, "Apply energy dispersion?"
  23 +refit, b, h, no,,, "Do refitting?"
  24 +caldb, s, a, "dummy",,, "Calibration database"
  25 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  26 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  27 +outmodel, f, a, "crab_results.xml",,, "Source model output file"
28 28 #
29 29 # Standard parameters
30 30 #====================
... ...
src/ctmodel/ctmodel.cpp
  View file @033ed45
... ... @@ -33,6 +33,7 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
  36 +#define G_GET_PARAMETERS "ctmodel::get_parameters()"
36 37  
37 38 /* __ Debug definitions __________________________________________________ */
38 39  
... ... @@ -207,9 +208,6 @@ void ctmodel::run(void)
207 208 // Get task parameters
208 209 get_parameters();
209 210  
210   - // Setup observation container
211   - setup_obs();
212   -
213 211 // Write parameters into logger
214 212 if (logTerse()) {
215 213 log_parameters();
... ... @@ -246,11 +244,6 @@ void ctmodel::run(void)
246 244 log.header1("Generate model cube");
247 245 }
248 246  
249   - // Initialise counts cube if cube has not yet been set
250   - if (!has_cube()) {
251   - init_cube();
252   - }
253   -
254 247 // Loop over all observations in the container
255 248 for (int i = 0; i < m_obs.size(); ++i) {
256 249  
... ... @@ -318,10 +311,10 @@ void ctmodel::save(void)
318 311 }
319 312  
320 313 // Make sure we have the FITS filename
321   - m_outfile = (*this)["outfile"].filename();
  314 + m_outcube = (*this)["outcube"].filename();
322 315  
323 316 // Save model cube into FITS file
324   - m_cube.save(m_outfile, clobber());
  317 + m_cube.save(m_outcube, clobber());
325 318  
326 319 // Return
327 320 return;
... ... @@ -365,35 +358,12 @@ void ctmodel::cube(const GCTAEventCube&amp; cube)
365 358 void ctmodel::init_members(void)
366 359 {
367 360 // Initialise members
368   - m_infile.clear();
369   - m_obsfile.clear();
370   - m_outfile.clear();
371   - m_caldb.clear();
372   - m_irf.clear();
373   - m_srcmdl.clear();
374   - m_proj.clear();
375   - m_coordsys.clear();
376   - m_ebinalg.clear();
377   - m_ebinfile.clear();
378   - m_ra = 0.0;
379   - m_dec = 0.0;
380   - m_deadc = 1.0;
381   - m_tmin = 0.0;
382   - m_tmax = 0.0;
383   - m_emin = 0.0;
384   - m_emax = 0.0;
385   - m_enumbins = 0;
386   - m_xref = 0.0;
387   - m_yref = 0.0;
388   - m_binsz = 0.0;
389   - m_nxpix = 0;
390   - m_nypix = 0;
  361 + m_outcube.clear();
391 362 m_apply_edisp = false;
392 363  
393 364 // Initialise protected members
394 365 m_obs.clear();
395 366 m_cube.clear();
396   - m_ebounds.clear();
397 367 m_gti.clear();
398 368 m_has_cube = false;
399 369  
... ... @@ -410,35 +380,12 @@ void ctmodel::init_members(void)
410 380 void ctmodel::copy_members(const ctmodel& app)
411 381 {
412 382 // Copy attributes
413   - m_infile = app.m_infile;
414   - m_obsfile = app.m_obsfile;
415   - m_outfile = app.m_outfile;
416   - m_caldb = app.m_caldb;
417   - m_irf = app.m_irf;
418   - m_srcmdl = app.m_srcmdl;
419   - m_proj = app.m_proj;
420   - m_coordsys = app.m_coordsys;
421   - m_ebinalg = app.m_ebinalg;
422   - m_ebinfile = app.m_ebinfile;
423   - m_ra = app.m_ra;
424   - m_dec = app.m_dec;
425   - m_deadc = app.m_deadc;
426   - m_tmin = app.m_tmin;
427   - m_tmax = app.m_tmax;
428   - m_emin = app.m_emin;
429   - m_emax = app.m_emax;
430   - m_enumbins = app.m_enumbins;
431   - m_xref = app.m_xref;
432   - m_yref = app.m_yref;
433   - m_binsz = app.m_binsz;
434   - m_nxpix = app.m_nxpix;
435   - m_nypix = app.m_nypix;
  383 + m_outcube = app.m_outcube;
436 384 m_apply_edisp = app.m_apply_edisp;
437 385  
438 386 // Copy protected members
439 387 m_obs = app.m_obs;
440 388 m_cube = app.m_cube;
441   - m_ebounds = app.m_ebounds;
442 389 m_gti = app.m_gti;
443 390 m_has_cube = app.m_has_cube;
444 391  
... ... @@ -465,255 +412,65 @@ void ctmodel::free_members(void)
465 412 ***************************************************************************/
466 413 void ctmodel::get_parameters(void)
467 414 {
468   - // Read input and optionally output cube filenames
469   - if (!has_cube()) {
470   - m_infile = (*this)["infile"].filename();
471   - }
472   - if (read_ahead()) {
473   - m_outfile = (*this)["outfile"].filename();
  415 + // If there are no observations in container then load them via user
  416 + // parameters
  417 + if (m_obs.size() == 0) {
  418 + m_obs = get_observations();
474 419 }
475 420  
476   - // Read model definiton file if required
  421 + // Read model definition file if required
477 422 if (m_obs.models().size() == 0) {
478   - m_srcmdl = (*this)["srcmdl"].filename();
  423 + std::string inmodel = (*this)["inmodel"].filename();
  424 +
  425 + // Load models from file
  426 + m_obs.models(inmodel);
  427 +
  428 + } // endif: there were no models
  429 +
  430 + // Read optionally output cube filenames
  431 + if (read_ahead()) {
  432 + m_outcube = (*this)["outcube"].filename();
479 433 }
480 434  
481 435 // Get energy dispersion flag parameters
482 436 m_apply_edisp = (*this)["edisp"].boolean();
483 437  
484   - // If the observation container is empty then read the parameters that
485   - // are necessary to fill it
486   - if (m_obs.size() == 0) {
487   -
488   - // Read observation definition filename and response parameters
489   - m_obsfile = (*this)["obsfile"].filename();
490   -
491   - // If no observation definition file has been specified then read all
492   - // parameters that are necessary to create an observation from scratch
493   - // (see method setup_obs)
494   - if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) {
495   - m_ra = (*this)["ra"].real();
496   - m_dec = (*this)["dec"].real();
497   - m_deadc = (*this)["deadc"].real();
498   - m_tmin = (*this)["tmin"].real();
499   - m_tmax = (*this)["tmax"].real();
500   - m_emin = (*this)["emin"].real();
501   - m_emax = (*this)["emax"].real();
502   - }
  438 + if (!m_has_cube) {
503 439  
504   - }
  440 + // Read cube definition file
  441 + std::string incube = (*this)["incube"].filename();
505 442  
506   - // Check if we need response information
507   - for (int i = 0; i < m_obs.size(); ++i) {
508   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
509   - if (obs != NULL) {
510   - if (!obs->has_response()) {
511   - m_caldb = (*this)["caldb"].string();
512   - m_irf = (*this)["irf"].string();
513   - break;
514   - }
515   - }
516   - }
  443 + // If no cube file has been specified then create a cube from
  444 + // the task parameters ...
  445 + if ((incube == "NONE") ||
  446 + (gammalib::strip_whitespace(incube) == "")) {
  447 +
  448 + // Create cube from scratch
  449 + m_cube = create_cube(m_obs);
517 450  
518   - // If no cube file has been specified then read all parameters that
519   - // are necessary to create the cube from scratch (see method init_cube)
520   - if ((!has_cube()) &&
521   - ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == ""))) {
522   - m_ebinalg = (*this)["ebinalg"].string();
523   - if (m_ebinalg == "FILE") {
524   - m_ebinfile = (*this)["ebinfile"].filename();
525 451 }
526   - else {
527   - m_emin = (*this)["emin"].real();
528   - m_emax = (*this)["emax"].real();
529   - m_enumbins = (*this)["enumbins"].integer();
530   - }
531   - m_proj = (*this)["proj"].string();
532   - m_coordsys = (*this)["coordsys"].string();
533   - m_xref = (*this)["xref"].real();
534   - m_yref = (*this)["yref"].real();
535   - m_binsz = (*this)["binsz"].real();
536   - m_nxpix = (*this)["nxpix"].integer();
537   - m_nypix = (*this)["nypix"].integer();
538   - }
539   -
540   - // Return
541   - return;
542   -}
543   -
544 452  
545   -/***********************************************************************//**
546   - * @brief Setup observation container
547   - *
548   - * This method makes sure that we have an observation container with an
549   - * associated model and that all CTA observations have a valid response
550   - * functions.
551   - *
552   - * If on input the observation container is empty, the method sets up an
553   - * observation container by either loading the information from an
554   - * observation definition XML file, by loading an event list or a counts
555   - * cube as a single observation, or by constructing an observation from
556   - * the scratch from the task parameters.
557   - *
558   - * If no model exists yet in the observation container, the method will load
559   - * the model from the file given by the task parameter "srcmdl".
560   - *
561   - * The method then loops over all observations in the container and makes
562   - * sure that all CTA observations have a valid response function. If no
563   - * response function exists so far, the response function will be set
564   - * using the "caldb" and "irf" task parameters.
565   - ***************************************************************************/
566   -void ctmodel::setup_obs(void)
567   -{
568   - // If there are no observations in the container then get them either
569   - // from the observation definition file or from the task parameters
570   - if (m_obs.size() == 0) {
571   -
572   - // If no observation definition file has been specified, then create
573   - // an observation container from the task parameters
574   - if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) {
575   -
576   - // Set pointing direction
577   - GCTAPointing pnt;
578   - GSkyDir skydir;
579   - skydir.radec_deg(m_ra, m_dec);
580   - pnt.dir(skydir);
581   -
582   - // Setup time interval covered by model
583   - GGti gti;
584   - GTime tmin(m_tmin);
585   - GTime tmax(m_tmax);
586   - gti.append(tmin, tmax);
587   -
588   - // Setup energy range covered by model
589   - GEnergy emin(m_emin, "TeV");
590   - GEnergy emax(m_emax, "TeV");
591   - GEbounds ebounds(1, emin, emax);
592   -
593   - // Allocate CTA observation and empty event list
594   - GCTAObservation obs;
595   - GCTAEventList list;
596   -
597   - // Set event list GTI and energy boundaries and append event list
598   - // to observation (we need in fact only the GTI for the
599   - // computations)
600   - list.gti(gti);
601   - list.ebounds(ebounds);
602   - obs.events(list);
603   -
604   - // Set CTA observation attributes
605   - obs.pointing(pnt);
606   - obs.ontime(gti.ontime());
607   - obs.livetime(gti.ontime()*m_deadc);
608   - obs.deadc(m_deadc);
609   -
610   - // Append CTA observation to container
611   - m_obs.append(obs);
612   -
613   - } // endif: created single CTA observation from task parameters
614   -
615   - // ... otherwise try to load information from the file
  453 + // ... otherwise load the cube from file and reset all bins
  454 + // to zero
616 455 else {
617 456  
618   - // First try to open the file as an event list or counts map
619   - try {
620   -
621   - // Allocate CTA observation
622   - GCTAObservation obs;
623   -
624   - // Load CTA observation from FITS file
625   - obs.load(m_obsfile);
626   -
627   - // Append CTA observation to container
628   - m_obs.append(obs);
629   -
630   - }
631   -
632   - // ... otherwise try to open as XML file
633   - catch (GException::fits_open_error &e) {
634   -
635   - // Load observations from XML file. This will throw
636   - // an exception if it fails.
637   - m_obs.load(m_obsfile);
  457 + // Load cube from given file
  458 + m_cube.load(incube);
638 459  
  460 + // Set all cube bins to zero
  461 + for (int i = 0; i < m_cube.size(); ++i) {
  462 + m_cube[i]->counts(0.0);
639 463 }
640   -
641   - } // endelse: loaded information from input file
642   -
643   - } // endif: there was no observation in the container
644   -
645   - // If there are no models associated with the observations then
646   - // load now the model definition from the XML file
647   - if (m_obs.models().size() == 0) {
648   - m_obs.models(GModels(m_srcmdl));
649   - }
650   -
651   - // For all observations that have no response, set the response
652   - // from the task parameters
653   - set_response(m_obs);
654   -
655   - // Return
656   - return;
657   -}
658   -
659   -
660   -/***********************************************************************//**
661   - * @brief Initialise model cube
662   - *
663   - * Initialises model cube by either reading the cube definition from an
664   - * existing cube file specified by the "infile" parameter, or by constructing
665   - * a cube from the energy and spatial binning information provided by the
666   - * task parameters. The GTIs of the model cube will be undefined on
667   - * exit.
668   - ***************************************************************************/
669   -void ctmodel::init_cube(void)
670   -{
671   - // Clear model cube and GTIs
672   - m_cube.clear();
673   - m_gti.clear();
674   -
675   - // If no cube file has been specified then create a cube from the task
676   - // parameters
677   - if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
678   -
679   - // Set dummy GTI that is needed for event cube creation
680   - GGti gti;
681   - gti.append(GTime(0.0), GTime(0.1234));
682   -
683   - // Set energy boundaries
684   - m_ebounds = get_ebounds();
685   -
686   - // Setup skymap
687   - GSkymap map = GSkymap(m_proj, m_coordsys,
688   - m_xref, m_yref, -m_binsz, m_binsz,
689   - m_nxpix, m_nypix, m_enumbins);
690   -
691   - // Create model cube
692   - m_cube = GCTAEventCube(map, m_ebounds, gti);
693   -
694   - }
695   -
696   - // ... otherwise load cube from file
697   - else {
698   -
699   - // Load cube
700   - m_cube.load(m_infile);
701   -
702   - // Set all cube bins to zero
703   - for (int i = 0; i < m_cube.size(); ++i) {
704   - m_cube[i]->counts(0.0);
705 464 }
706 465  
  466 + // Signal that cube has been set
  467 + m_has_cube = true;
707 468 }
708 469  
709   - // Signal that cube has been set
710   - m_has_cube = true;
711   -
712 470 // Return
713 471 return;
714 472 }
715 473  
716   -
717 474 /***********************************************************************//**
718 475 * @brief Fill model into model cube
719 476 *
... ...
src/ctmodel/ctmodel.hpp
  View file @033ed45
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTMODEL_NAME "ctmodel"
37   -#define CTMODEL_VERSION "00-03-00"
  37 +#define CTMODEL_VERSION "00-04-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -78,41 +78,16 @@ protected:
78 78 void copy_members(const ctmodel& app);
79 79 void free_members(void);
80 80 void get_parameters(void);
81   - void setup_obs(void);
82   - void init_cube(void);
83 81 void fill_cube(const GCTAObservation* obs);
84 82 bool has_cube(void) const;
85 83  
86 84 // User parameters
87   - std::string m_infile; //!< Input counts cube
88   - std::string m_obsfile; //!< Event list, counts map or observation definition XML file
89   - std::string m_outfile; //!< Output model cube
90   - std::string m_caldb; //!< Calibration database
91   - std::string m_irf; //!< Instrument response functions
92   - std::string m_srcmdl; //!< Source model
93   - double m_ra; //!< RA of pointing direction
94   - double m_dec; //!< DEC of pointing direction
95   - double m_deadc; //!< Deadtime correction factor
96   - double m_tmin; //!< Start time
97   - double m_tmax; //!< Stop time
98   - std::string m_ebinalg; //!< Algorithm for energy binning
99   - std::string m_ebinfile; //!< FITS-file containing energy binning
100   - double m_emin; //!< Lower energy
101   - double m_emax; //!< Upper energy
102   - int m_enumbins; //!< Number of energy bins
103   - std::string m_proj; //!< WCS projection
104   - std::string m_coordsys; //!< Coordinate system
105   - double m_xref; //!< Longitude reference coordinate
106   - double m_yref; //!< Latitude reference coordinate
107   - double m_binsz; //!< Pixel size
108   - int m_nxpix; //!< Number of pixels in longitude
109   - int m_nypix; //!< Number of pixels in latitude
  85 + std::string m_outcube; //!< Output model cube
110 86 bool m_apply_edisp; //!< Apply energy dispersion?
111 87  
112 88 // Protected members
113 89 GObservations m_obs; //!< Observation container
114 90 GCTAEventCube m_cube; //!< Model cube
115   - GEbounds m_ebounds; //!< Energy boundaries
116 91 GGti m_gti; //!< Model cube GTIs
117 92 bool m_has_cube; //!< Signal if cube has been set or loaded
118 93 };
... ...
src/ctmodel/ctmodel.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctmodel.par #
8   -# Version: 00-03-00 #
  8 +# Version: 00-04-00 #
9 9 # #
10 10 # Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -17,18 +17,19 @@
17 17 #
18 18 # General parameters
19 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   -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  20 +incube, f, a, "NONE",,, "Input counts map or observation definition file"
  21 +outcube, f, a, "modmap.fits",,, "Output counts map or observation definition file"
  22 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
23 23 edisp, b, h, no,,, "Apply energy dispersion?"
24 24 #
25 25 # Observation parameters
26 26 #=======================
27   -obsfile, f, a, "NONE",,, "Observation definition file"
  27 +inobs, f, a, "NONE",,, "Input event list, cube or observation definition file"
28 28 caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
29 29 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
30 30 ra, r, a, 83.63,0,360, "RA of pointing (degrees)"
31 31 dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)"
  32 +rad, r, a, 5.0, 0,180, "Radius of FOV (degrees)"
32 33 tmin, r, a, 0.0,0,, "Start time (MET in s)"
33 34 tmax, r, a, 1800.0,0,, "End time (MET in s)"
34 35 deadc, r, h, 0.95,0,1, "Average deadtime correction factor"
... ... @@ -43,6 +44,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
43 44 #
44 45 # Spatial binning parameters
45 46 #===========================
  47 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
46 48 nxpix, i, a, 200,,, "Size of the X axis in pixels"
47 49 nypix, i, a, 200,,, "Size of the Y axis in pixels"
48 50 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/ctobssim/ctobssim.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctobssim - Observation simulator tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -377,7 +377,7 @@ void ctobssim::run(void)
377 377 outfile = m_prefix + gammalib::str(i) + ".fits";
378 378 }
379 379 else {
380   - outfile = (*this)["outfile"].filename();
  380 + outfile = (*this)["outevents"].filename();
381 381 }
382 382  
383 383 // Store output file name in original observation
... ... @@ -490,35 +490,20 @@ void ctobssim::save(void)
490 490 void ctobssim::init_members(void)
491 491 {
492 492 // Initialise user parameters
493   - m_infile.clear();
494   - m_outfile.clear();
  493 + m_outevents.clear();
495 494 m_prefix.clear();
496   - m_seed = 1;
497   - m_ra = 0.0;
498   - m_dec = 0.0;
499   - m_rad = 0.0;
500   - m_tmin = 0.0;
501   - m_tmax = 0.0;
502   - m_emin = 0.0;
503   - m_emax = 0.0;
  495 + m_seed = 1;
504 496 m_apply_edisp = false;
505   - m_deadc = 1.0;
506 497  
507 498 // Initialise protected members
508 499 m_rans.clear();
509 500 m_obs.clear();
510   - m_use_xml = false;
511 501 m_save_and_dispose = false;
512 502  
513 503 // Set fixed parameters
514 504 m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2
515 505 m_max_photons = 1000000; //!< Maximum number of photons / time slice
516 506  
517   - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
518   - // which is defined in GCTALib.hpp. This is somehow a kluge. We need
519   - // a better mechanism to implement the CTA reference MJD.
520   - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
521   -
522 507 // Initialise first event identifier
523 508 m_event_id = 1;
524 509  
... ... @@ -535,28 +520,17 @@ void ctobssim::init_members(void)
535 520 void ctobssim::copy_members(const ctobssim& app)
536 521 {
537 522 // Copy user parameters
538   - m_infile = app.m_infile;
539   - m_outfile = app.m_outfile;
  523 + m_outevents = app.m_outevents;
540 524 m_prefix = app.m_prefix;
541 525 m_seed = app.m_seed;
542   - m_ra = app.m_ra;
543   - m_dec = app.m_dec;
544   - m_rad = app.m_rad;
545   - m_tmin = app.m_tmin;
546   - m_tmax = app.m_tmax;
547   - m_emin = app.m_emin;
548   - m_emax = app.m_emax;
549 526 m_apply_edisp = app.m_apply_edisp;
550   - m_deadc = app.m_deadc;
551 527  
552 528 // Copy protected members
553 529 m_area = app.m_area;
554 530 m_max_photons = app.m_max_photons;
555 531 m_rans = app.m_rans;
556 532 m_obs = app.m_obs;
557   - m_use_xml = app.m_use_xml;
558 533 m_save_and_dispose = app.m_save_and_dispose;
559   - m_cta_ref = app.m_cta_ref;
560 534 m_event_id = app.m_event_id;
561 535  
562 536 // Return
... ... @@ -585,47 +559,16 @@ void ctobssim::free_members(void)
585 559 ***************************************************************************/
586 560 void ctobssim::get_parameters(void)
587 561 {
588   - // If there are no observations in container then add a single CTA
589   - // observation using the parameters from the parameter file
  562 + // If there are no observations in container then load them via user
  563 + // parameters
590 564 if (m_obs.size() == 0) {
  565 + m_obs = get_observations();
  566 + }
591 567  
592   - // Get CTA observation parameters
593   - m_infile = (*this)["infile"].filename();
594   - m_ra = (*this)["ra"].real();
595   - m_dec = (*this)["dec"].real();
596   -
597   - // Set pointing direction
598   - GCTAPointing pnt;
599   - GSkyDir skydir;
600   - skydir.radec_deg(m_ra, m_dec);
601   - pnt.dir(skydir);
602   -
603   - // Allocate CTA observation
604   - GCTAObservation obs;
605   -
606   - // Set response
607   - set_obs_response(&obs);
608   -
609   - // Set pointing
610   - obs.pointing(pnt);
611   -
612   - // Set event list (queries remaining parameters)
613   - set_list(&obs);
614   -
615   - // Append CTA observation to container
616   - m_obs.append(obs);
617   -
618   - // Load models into container
619   - m_obs.models(m_infile);
620   -
621   - // Signal that no XML file should be used for storage
622   - m_use_xml = false;
623   -
624   - } // endif: there was no observation in the container
625   -
626   - // ... otherwise we signal that an XML file should be written
627   - else {
628   - m_use_xml = true;
  568 + // Read model definition file if required
  569 + if (m_obs.models().size() == 0) {
  570 + std::string inmodel = (*this)["inmodel"].filename();
  571 + m_obs.models(inmodel);
629 572 }
630 573  
631 574 // Get other parameters
... ... @@ -635,18 +578,8 @@ void ctobssim::get_parameters(void)
635 578 // Optionally read ahead parameters so that they get correctly
636 579 // dumped into the log file
637 580 if (read_ahead()) {
638   - m_outfile = (*this)["outfile"].filename();
639   - m_prefix = (*this)["prefix"].string();
640   - }
641   -
642   - // Make sure that all observations hold a CTA event list
643   - for (int i = 0; i < m_obs.size(); ++i) {
644   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
645   - if (obs != NULL) {
646   - if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
647   - set_list(obs);
648   - }
649   - }
  581 + m_outevents = (*this)["outevents"].filename();
  582 + m_prefix = (*this)["prefix"].string();
650 583 }
651 584  
652 585 // Initialise random number generators. We initialise here one random
... ... @@ -700,81 +633,6 @@ void ctobssim::get_parameters(void)
700 633  
701 634  
702 635 /***********************************************************************//**
703   - * @brief Set empty CTA event list
704   - *
705   - * @param[in] obs CTA observation.
706   - *
707   - * Attaches an empty event list to CTA observation. The method also sets the
708   - * pointing direction using the m_ra and m_dec members, the ROI based on
709   - * m_ra, m_dec and m_rad, a single GTI based on m_tmin and m_tmax, and a
710   - * single energy boundary based on m_emin and m_emax. The method furthermore
711   - * sets the ontime, livetime and deadtime correction factor.
712   - ***************************************************************************/
713   -void ctobssim::set_list(GCTAObservation* obs)
714   -{
715   - // Continue only if observation is valid
716   - if (obs != NULL) {
717   -
718   - // Get CTA observation parameters
719   - m_ra = (*this)["ra"].real();
720   - m_dec = (*this)["dec"].real();
721   - m_rad = (*this)["rad"].real();
722   - m_tmin = (*this)["tmin"].real();
723   - m_tmax = (*this)["tmax"].real();
724   - m_emin = (*this)["emin"].real();
725   - m_emax = (*this)["emax"].real();
726   - m_deadc = (*this)["deadc"].real();
727   -
728   - // Allocate CTA event list
729   - GCTAEventList events;
730   -
731   - // Set pointing direction
732   - GCTAPointing pnt;
733   - GSkyDir skydir;
734   - skydir.radec_deg(m_ra, m_dec);
735   - pnt.dir(skydir);
736   -
737   - // Set ROI
738   - GCTAInstDir instdir(skydir);
739   - GCTARoi roi(instdir, m_rad);
740   -
741   - // Set GTI
742   - GGti gti(m_cta_ref);
743   - GTime tstart;
744   - GTime tstop;
745   - tstart.set(m_tmin, m_cta_ref);
746   - tstop.set(m_tmax, m_cta_ref);
747   - gti.append(tstart, tstop);
748   -
749   - // Set energy boundaries
750   - GEbounds ebounds;
751   - GEnergy emin;
752   - GEnergy emax;
753   - emin.TeV(m_emin);
754   - emax.TeV(m_emax);
755   - ebounds.append(emin, emax);
756   -
757   - // Set CTA event list attributes
758   - events.roi(roi);
759   - events.gti(gti);
760   - events.ebounds(ebounds);
761   -
762   - // Attach event list to CTA observation
763   - obs->events(events);
764   -
765   - // Set observation ontime, livetime and deadtime correction factor
766   - obs->ontime(gti.ontime());
767   - obs->livetime(gti.ontime()*m_deadc);
768   - obs->deadc(m_deadc);
769   -
770   - } // endif: oberservation was valid
771   -
772   - // Return
773   - return;
774   -}
775   -
776   -
777   -/***********************************************************************//**
778 636 * @brief Simulate source events from photon list
779 637 *
780 638 * @param[in] obs Pointer on CTA observation.
... ... @@ -1163,13 +1021,13 @@ void ctobssim::save_fits(void)
1163 1021 if (!m_save_and_dispose) {
1164 1022  
1165 1023 // Get output filename
1166   - m_outfile = (*this)["outfile"].filename();
  1024 + m_outevents = (*this)["outevents"].filename();
1167 1025  
1168 1026 // Get CTA observation from observation container
1169 1027 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1170 1028  
1171 1029 // Save observation into FITS file
1172   - obs->save(m_outfile, clobber());
  1030 + obs->save(m_outevents, clobber());
1173 1031  
1174 1032 } // endif: event list has not yet been saved and disposed
1175 1033  
... ... @@ -1191,14 +1049,14 @@ void ctobssim::save_fits(void)
1191 1049 void ctobssim::save_xml(void)
1192 1050 {
1193 1051 // Get output filename and prefix
1194   - m_outfile = (*this)["outfile"].filename();
  1052 + m_outevents = (*this)["outevents"].filename();
1195 1053 m_prefix = (*this)["prefix"].string();
1196 1054  
1197 1055 // Issue warning if output filename has no .xml suffix
1198   - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
  1056 + std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4));
1199 1057 if (suffix != ".xml") {
1200 1058 log << "*** WARNING: Name of observation definition output file \""+
1201   - m_outfile+"\"" << std::endl;
  1059 + m_outevents+"\"" << std::endl;
1202 1060 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1203 1061 log << "*** WARNING: This is not an error, but might be misleading."
1204 1062 " It is recommended" << std::endl;
... ... @@ -1239,7 +1097,7 @@ void ctobssim::save_xml(void)
1239 1097 } // endif: event list has not yet been saved and disposed
1240 1098  
1241 1099 // Save observations in XML file
1242   - m_obs.save(m_outfile);
  1100 + m_obs.save(m_outevents);
1243 1101  
1244 1102 // Return
1245 1103 return;
... ...
src/ctobssim/ctobssim.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctobssim - Observation simulator tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTOBSSIM_NAME "ctobssim"
37   -#define CTOBSSIM_VERSION "00-09-00"
  37 +#define CTOBSSIM_VERSION "00-10-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -73,7 +73,6 @@ protected:
73 73 void copy_members(const ctobssim& app);
74 74 void free_members(void);
75 75 void get_parameters(void);
76   - void set_list(GCTAObservation* obs);
77 76 void simulate_source(GCTAObservation* obs,
78 77 const GModels& models,
79 78 GRan& ran,
... ... @@ -86,19 +85,10 @@ protected:
86 85 void save_xml(void);
87 86  
88 87 // User parameters
89   - std::string m_infile; //!< Input model
90   - std::string m_outfile; //!< Output events file
  88 + std::string m_outevents; //!< Output events file
91 89 std::string m_prefix; //!< Prefix for multiple event lists
92 90 int m_seed; //!< Random number generator seed
93   - double m_ra; //!< RA of pointing direction
94   - double m_dec; //!< DEC of pointing direction
95   - double m_rad; //!< FOV radius (degrees)
96   - double m_tmin; //!< Start time (MET)
97   - double m_tmax; //!< Stop time (MET)
98   - double m_emin; //!< Lower energy (TeV)
99   - double m_emax; //!< Upper energy (TeV)
100 91 bool m_apply_edisp; //!< Apply energy dispersion?
101   - double m_deadc; //!< Average deadtime correction
102 92  
103 93 // Protected members
104 94 mutable GObservations m_obs; //!< Observation container
... ... @@ -106,8 +96,6 @@ protected:
106 96 double m_area; //!< Surface area for simulation (cm2)
107 97 int m_max_photons; //!< Maximum number of photons/slice
108 98 std::vector<GRan> m_rans; //!< Random number generators
109   - bool m_use_xml; //!< Use XML file instead of FITS file
110   - GTimeReference m_cta_ref; //!< CTA time reference
111 99 int m_event_id; //!< Event identifier
112 100 };
113 101  
... ... @@ -120,7 +108,6 @@ protected:
120 108 inline
121 109 const GObservations& ctobssim::obs(void) const
122 110 {
123   - // Return observation container
124 111 return m_obs;
125 112 }
126 113  
... ...
src/ctobssim/ctobssim.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctobssim.par #
8   -# Version: 00-08-00 #
  8 +# Version: 00-10-00 #
9 9 # #
10 10 # Author: Jurgen Knodlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -17,9 +17,10 @@
17 17 #
18 18 # File information
19 19 #=================
20   -infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model"
21   -outfile,f, a, "events.fits",,, "Output event data file or observation definition file"
22   -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file"
  20 +inobs, f, h, "NONE",,,"Input event list or observation definition file"
  21 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Model"
  22 +outevents, f, a, "events.fits",,,"Output event data file or observation definition file"
  23 +prefix, s, h, "sim_events_",,,"Prefix for event lists in observation definition file"
23 24 #
24 25 # Observation parameters
25 26 #=======================
... ...
src/ctpsfcube/ctpsfcube.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctpsfcube - PSF cube generation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Chia-Chun Lu *
  4 + * copyright (C) 2014-2015 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,7 +35,6 @@
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
37 37 #define G_GET_PARAMETERS "ctpsfcube::get_parameters()"
38   -#define G_SET_FROM_CNTMAP "ctpsfcube::set_from_cntmap(std::string&)"
39 38  
40 39 /* __ Debug definitions __________________________________________________ */
41 40  
... ... @@ -275,10 +274,10 @@ void ctpsfcube::save(void)
275 274 }
276 275  
277 276 // Get output filename
278   - m_outfile = (*this)["outfile"].filename();
  277 + m_outcube = (*this)["outcube"].filename();
279 278  
280 279 // Save PSF cube
281   - m_psfcube.save(m_outfile, clobber());
  280 + m_psfcube.save(m_outcube, clobber());
282 281  
283 282 // Return
284 283 return;
... ... @@ -297,13 +296,12 @@ void ctpsfcube::save(void)
297 296 void ctpsfcube::init_members(void)
298 297 {
299 298 // Initialise members
300   - m_outfile.clear();
  299 + m_outcube.clear();
301 300 m_apply_edisp = false;
302 301  
303 302 // Initialise protected members
304 303 m_obs.clear();
305 304 m_psfcube.clear();
306   - m_ebounds.clear();
307 305  
308 306 // Return
309 307 return;
... ... @@ -318,13 +316,12 @@ void ctpsfcube::init_members(void)
318 316 void ctpsfcube::copy_members(const ctpsfcube& app)
319 317 {
320 318 // Copy attributes
321   - m_outfile = app.m_outfile;
  319 + m_outcube = app.m_outcube;
322 320 m_apply_edisp = app.m_apply_edisp;
323 321  
324 322 // Copy protected members
325 323 m_obs = app.m_obs;
326 324 m_psfcube = app.m_psfcube;
327   - m_ebounds = app.m_ebounds;
328 325  
329 326 // Return
330 327 return;
... ... @@ -346,167 +343,56 @@ void ctpsfcube::free_members(void)
346 343 *
347 344 * Get all task parameters from parameter file or (if required) by querying
348 345 * the user. The parameters are read in the correct order.
349   - *
350   - * @todo Setup PSF cube from counts map
351 346 ***************************************************************************/
352 347 void ctpsfcube::get_parameters(void)
353 348 {
354   - // If we do not have any observations in the container then get an
355   - // input file name or observation descriptor file
  349 + // If there are no observations in container then load them via user
  350 + // parameters
356 351 if (m_obs.size() == 0) {
357   - get_obs();
358   - }
359   -
360   - // Make sure that response is set
361   - set_response(m_obs);
362 352  
363   - // Read energy dispersion flag
364   - m_apply_edisp = (*this)["edisp"].boolean();
  353 + // Throw exception if no input observation file is given
  354 + require_inobs(G_GET_PARAMETERS);
365 355  
366   - // If no counts map is specified then setup the PSF cube from
367   - // the user parameters
368   - std::string cntmap = (*this)["cntmap"].filename();
369   - if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) {
370   -
371   - // Get user parameters for counts map definition
372   - std::string wcs = (*this)["proj"].string();
373   - std::string coordsys = (*this)["coordsys"].string();
374   - double xref = (*this)["xref"].real();
375   - double yref = (*this)["yref"].real();
376   - double binsz = (*this)["binsz"].real();
377   - int nxpix = (*this)["nxpix"].integer();
378   - int nypix = (*this)["nypix"].integer();
379   - double dmax = (*this)["amax"].real();
380   - int ndbins = (*this)["anumbins"].integer();
381   -
382   - // Get energy definition
383   - m_ebounds = get_ebounds();
384   -
385   - // Define PSF cube
386   - m_psfcube = GCTAMeanPsf(wcs, coordsys, xref, yref,
387   - -binsz, binsz, nxpix, nypix,
388   - m_ebounds, dmax, ndbins);
389   -
390   - } // endif: PSF cube set from user parameters
391   -
392   - // ... otherwise setup the PSF cube from the counts map
393   - else {
394   -
395   - // Set PSF cube from counts map
396   - set_from_cntmap(cntmap);
397   -
398   - }
399   -
400   - // Read output filename (if needed)
401   - if (read_ahead()) {
402   - m_outfile = (*this)["outfile"].filename();
403   - }
  356 + // Build observation container
  357 + m_obs = get_observations();
404 358  
405   - // Return
406   - return;
407   -}
  359 + } // endif: there was no observation in the container
408 360  
  361 + // Get the incube filename
  362 + std::string incube = (*this)["incube"].filename();
409 363  
410   -/***********************************************************************//**
411   - * @brief Get observation definition
412   - *
413   - * Get observation definition from the user parameters.
414   - ***************************************************************************/
415   -void ctpsfcube::get_obs(void)
416   -{
417   - // Get input filename
418   - std::string filename = (*this)["infile"].filename();
  364 + // Get additional binning parameters
  365 + double amax = (*this)["amax"].real();
  366 + int anumbins = (*this)["anumbins"].integer();
419 367  
420   - // Try first to open as FITS file
421   - try {
  368 + // Check for filename validity
  369 + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
422 370  
423   - // Allocate CTA observation
424   - GCTAObservation obs;
425   -
426   - // Load input file in CTA observation
427   - obs.load(filename);
  371 + // Create an event cube based on task parameters
  372 + GCTAEventCube cube = create_cube(m_obs);
428 373  
429   - // Append CTA observation to container
430   - m_obs.append(obs);
  374 + // Define psf cube
  375 + m_psfcube = GCTAMeanPsf(cube, amax, anumbins);
431 376  
432 377 }
433   -
434   - // ... otherwise try to open as XML file
435   - catch (GException::fits_open_error &e) {
436 378  
437   - // Load observations from XML file
438   - m_obs.load(filename);
  379 + // ... otherwise setup the exposure cube from the counts map
  380 + else {
439 381  
440   - }
  382 + // Load event cube from filename
  383 + GCTAEventCube cube(incube);
441 384  
442   - // Return
443   - return;
444   -}
  385 + // Define psf cube
  386 + m_psfcube = GCTAMeanPsf(cube, amax, anumbins);
445 387  
  388 + } // endelse: cube loaded from file
446 389  
447   -/***********************************************************************//**
448   - * @brief Set PSF cube definition from counts map
449   - *
450   - * @exception GException::invalid_value
451   - * Invalid counts map projection or invalid events encountered.
452   - *
453   - * Set PSF cube definition from counts map.
454   - ***************************************************************************/
455   -void ctpsfcube::set_from_cntmap(const std::string& filename)
456   -{
457   - // Allocate CTA observation
458   - GCTAObservation obs;
459   -
460   - // Load counts map in CTA observation
461   - obs.load(filename);
462   -
463   - // Set PSF cube from counts map
464   - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
465   -
466   - // Continue only if cube is valid
467   - if (cube != NULL) {
468   -
469   - // Get sky map projection
470   - const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
471   -
472   - // Continue only if projection is valid
473   - if (wcs != NULL) {
474   -
475   - // Get user parameters for counts map definition
476   - std::string proj = wcs->code();
477   - std::string coordsys = wcs->coordsys();
478   - double xref = wcs->crval(0);
479   - double yref = wcs->crval(1);
480   - double dx = wcs->cdelt(0);
481   - double dy = wcs->cdelt(1);
482   - int nx = cube->map().nx();
483   - int ny = cube->map().ny();
484   - double amax = (*this)["amax"].real();
485   - int anumbins = (*this)["anumbins"].integer();
486   -
487   - // Get energy definition
488   - m_ebounds = cube->ebounds();
489   -
490   - // Define PSF cube
491   - m_psfcube = GCTAMeanPsf(proj, coordsys, xref, yref,
492   - dx, dy, nx, ny,
493   - m_ebounds, amax, anumbins);
494   -
495   - } // endif: WCS projection was valid
496   -
497   - // ... projection is not of WCS type
498   - else {
499   - std::string msg = "Counts map project is not of WCS type.";
500   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
501   - }
502   -
503   - } // endif: observation contained an events cube
  390 + // Read energy dispersion flag
  391 + m_apply_edisp = (*this)["edisp"].boolean();
504 392  
505   - // ... there is not events cube
506   - else {
507   - std::string msg = "No events cube found in file \""
508   - ""+filename+"\".";
509   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  393 + // Read output filename (if needed)
  394 + if (read_ahead()) {
  395 + m_outcube = (*this)["outcube"].filename();
510 396 }
511 397  
512 398 // Return
... ...
src/ctpsfcube/ctpsfcube.hpp
  View file @033ed45
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTPSFCUBE_NAME "ctpsfcube"
37   -#define CTPSFCUBE_VERSION "00-01-00"
  37 +#define CTPSFCUBE_VERSION "00-02-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -67,17 +67,14 @@ protected:
67 67 void copy_members(const ctpsfcube& app);
68 68 void free_members(void);
69 69 void get_parameters(void);
70   - void get_obs(void);
71   - void set_from_cntmap(const std::string& filename);
72 70  
73 71 // User parameters
74   - std::string m_outfile; //!< Output PSF cube file
  72 + std::string m_outcube; //!< Output PSF cube file
75 73 bool m_apply_edisp; //!< Apply energy dispersion?
76 74  
77 75 // Protected members
78 76 GObservations m_obs; //!< Observation container
79 77 GCTAMeanPsf m_psfcube; //!< PSF cube
80   - GEbounds m_ebounds; //!< Energy boundaries
81 78 };
82 79  
83 80  
... ...
src/ctpsfcube/ctpsfcube.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctpsfcube.par #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Chia-Chun Lu #
11 11 # c.lu@mpi-hd.mpg.de #
... ... @@ -15,11 +15,11 @@
15 15 # #
16 16 #############################################################################
17 17 #
18   -# General parameters
  18 +# Observational parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Event list or observation definition file"
21   -cntmap, f, a, "NONE",,, "Counts map for exposure cube definition"
22   -outfile, f, a, "psfcube.fits",,, "Output exposure cube file"
  20 +inobs, f, a, "NONE",,, "Event list or observation definition file"
  21 +incube, f, a, "NONE",,, "Counts map for psf cube definition"
  22 +outcube, f, a, "psfcube.fits",,, "Output psf cube file"
23 23 caldb, s, a, "dummy",,, "Calibration database"
24 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 25 edisp, b, h, no,,, "Apply energy dispersion?"
... ... @@ -34,6 +34,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
34 34 #
35 35 # Spatial binning parameters
36 36 #===========================
  37 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
37 38 nxpix, i, a, 200,,, "Size of the X axis in pixels"
38 39 nypix, i, a, 200,,, "Size of the Y axis in pixels"
39 40 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/ctselect/ctselect.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctselect - Data selection tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -387,8 +387,7 @@ void ctselect::save(void)
387 387 void ctselect::init_members(void)
388 388 {
389 389 // Initialise parameters
390   - m_infile.clear();
391   - m_outfile.clear();
  390 + m_outobs.clear();
392 391 m_prefix.clear();
393 392 m_usepnt = false;
394 393 m_ra = -1.0;
... ... @@ -406,12 +405,6 @@ void ctselect::init_members(void)
406 405 m_infiles.clear();
407 406 m_timemin.clear();
408 407 m_timemax.clear();
409   - m_use_xml = false;
410   -
411   - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
412   - // which is defined in GCTALib.hpp. This is somehow a kluge. We need
413   - // a better mechanism to implement the CTA reference MJD.
414   - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
415 408  
416 409 // Return
417 410 return;
... ... @@ -426,8 +419,7 @@ void ctselect::init_members(void)
426 419 void ctselect::copy_members(const ctselect& app)
427 420 {
428 421 // Copy parameters
429   - m_infile = app.m_infile;
430   - m_outfile = app.m_outfile;
  422 + m_outobs = app.m_outobs;
431 423 m_prefix = app.m_prefix;
432 424 m_usepnt = app.m_usepnt;
433 425 m_ra = app.m_ra;
... ... @@ -445,8 +437,6 @@ void ctselect::copy_members(const ctselect&amp; app)
445 437 m_infiles = app.m_infiles;
446 438 m_timemin = app.m_timemin;
447 439 m_timemax = app.m_timemax;
448   - m_use_xml = app.m_use_xml;
449   - m_cta_ref = app.m_cta_ref;
450 440  
451 441 // Return
452 442 return;
... ... @@ -466,7 +456,7 @@ void ctselect::free_members(void)
466 456 /***********************************************************************//**
467 457 * @brief Get application parameters
468 458 *
469   - * Get all task parameters from parameter file or (if required) by querying
  459 + * Get all user parameters from parameter file or (if required) by querying
470 460 * the user. Times are assumed to be in the native CTA MJD format.
471 461 *
472 462 * This method also loads observations if no observations are yet allocated.
... ... @@ -475,40 +465,15 @@ void ctselect::free_members(void)
475 465 ***************************************************************************/
476 466 void ctselect::get_parameters(void)
477 467 {
478   - // If there are no observations in container then add a single CTA
479   - // observation using the parameters from the parameter file
  468 + // If there are no observations in container then load them via user
  469 + // parameters
480 470 if (m_obs.size() == 0) {
481 471  
482   - // Get CTA event list file name
483   - m_infile = (*this)["infile"].filename();
484   -
485   - // Allocate CTA observation
486   - GCTAObservation obs;
487   -
488   - // Try first to open as FITS file
489   - try {
490   -
491   - // Load event list in CTA observation
492   - obs.load(m_infile);
493   -
494   - // Append CTA observation to container
495   - m_obs.append(obs);
496   -
497   - // Signal that no XML file should be used for storage
498   - m_use_xml = false;
499   -
500   - }
501   -
502   - // ... otherwise try to open as XML file
503   - catch (GException::fits_open_error &e) {
  472 + // Throw exception if no input observation file is given
  473 + require_inobs(G_GET_PARAMETERS);
504 474  
505   - // Load observations from XML file
506   - m_obs.load(m_infile);
507   -
508   - // Signal that XML file should be used for storage
509   - m_use_xml = true;
510   -
511   - }
  475 + // Get observation container without response (not needed)
  476 + m_obs = get_observations(false);
512 477  
513 478 } // endif: there was no observation in the container
514 479  
... ... @@ -537,8 +502,8 @@ void ctselect::get_parameters(void)
537 502 // Optionally read ahead parameters so that they get correctly
538 503 // dumped into the log file
539 504 if (read_ahead()) {
540   - m_outfile = (*this)["outfile"].filename();
541   - m_prefix = (*this)["prefix"].string();
  505 + m_outobs = (*this)["outobs"].filename();
  506 + m_prefix = (*this)["prefix"].string();
542 507 }
543 508  
544 509 // Set time interval with input times given in CTA reference
... ... @@ -1048,7 +1013,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const
1048 1013 if (missing.size() > 0) {
1049 1014 message = "The following columns are missing in the"
1050 1015 " \"EVENTS\" extension of input file \""
1051   - + m_outfile + "\": ";
  1016 + + m_outobs + "\": ";
1052 1017 for (int i = 0; i < missing.size(); ++i) {
1053 1018 message += "\"" + missing[i] + "\"";
1054 1019 if (i < missing.size()-1) {
... ... @@ -1060,7 +1025,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const
1060 1025 }
1061 1026 catch (GException::fits_hdu_not_found& e) {
1062 1027 message = "No \"EVENTS\" extension found in input file \""
1063   - + m_outfile + "\".";
  1028 + + m_outobs + "\".";
1064 1029 }
1065 1030  
1066 1031 // Return
... ... @@ -1103,13 +1068,13 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const
1103 1068 void ctselect::save_fits(void)
1104 1069 {
1105 1070 // Get output filename
1106   - m_outfile = (*this)["outfile"].filename();
  1071 + m_outobs = (*this)["outobs"].filename();
1107 1072  
1108 1073 // Get CTA observation from observation container
1109 1074 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1110 1075  
1111 1076 // Save event list
1112   - save_event_list(obs, m_infiles[0], m_outfile);
  1077 + save_event_list(obs, m_infiles[0], m_outobs);
1113 1078  
1114 1079 // Return
1115 1080 return;
... ... @@ -1129,14 +1094,14 @@ void ctselect::save_fits(void)
1129 1094 void ctselect::save_xml(void)
1130 1095 {
1131 1096 // Get output filename and prefix
1132   - m_outfile = (*this)["outfile"].filename();
  1097 + m_outobs = (*this)["outobs"].filename();
1133 1098 m_prefix = (*this)["prefix"].string();
1134 1099  
1135 1100 // Issue warning if output filename has no .xml suffix
1136   - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
  1101 + std::string suffix = gammalib::tolower(m_outobs.substr(m_outobs.length()-4,4));
1137 1102 if (suffix != ".xml") {
1138 1103 log << "*** WARNING: Name of observation definition output file \""+
1139   - m_outfile+"\"" << std::endl;
  1104 + m_outobs+"\"" << std::endl;
1140 1105 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1141 1106 log << "*** WARNING: This is not an error, but might be misleading."
1142 1107 " It is recommended" << std::endl;
... ... @@ -1167,7 +1132,7 @@ void ctselect::save_xml(void)
1167 1132 } // endfor: looped over observations
1168 1133  
1169 1134 // Save observations in XML file
1170   - m_obs.save(m_outfile);
  1135 + m_obs.save(m_outobs);
1171 1136  
1172 1137 // Return
1173 1138 return;
... ...
src/ctselect/ctselect.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctselect - Data selection tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,7 +36,7 @@
36 36  
37 37 /* __Definitions _________________________________________________________ */
38 38 #define CTSELECT_NAME "ctselect"
39   -#define CTSELECT_VERSION "00-06-00"
  39 +#define CTSELECT_VERSION "00-07-00"
40 40  
41 41  
42 42 /***********************************************************************//**
... ... @@ -82,8 +82,7 @@ protected:
82 82 const std::string& outfile) const;
83 83  
84 84 // User parameters
85   - std::string m_infile; //!< Input event list or XML file
86   - std::string m_outfile; //!< Output event list or XML file
  85 + std::string m_outobs; //!< Output event list or XML file
87 86 std::string m_prefix; //!< Prefix for multiple event lists
88 87 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
89 88 double m_ra; //!< RA of ROI centre
... ... @@ -101,8 +100,7 @@ protected:
101 100 std::vector<std::string> m_infiles; //!< Input event filenames
102 101 GTime m_timemin; //!< Earliest time
103 102 GTime m_timemax; //!< Latest time
104   - bool m_use_xml; //!< Use XML file instead of FITS file
105   - GTimeReference m_cta_ref; //!< CTA time reference
  103 +
106 104 };
107 105  
108 106  
... ...
src/ctselect/ctselect.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctselect.par #
8   -# Version: 00-05-00 #
  8 +# Version: 00-06-00 #
9 9 # #
10 10 # Author: Jurgen Knodlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -17,9 +17,11 @@
17 17 #
18 18 # File information
19 19 #=================
20   -infile, f, a, "events.fits",,,"Input event list or observation definition file"
21   -outfile,f, a, "selected_events.fits",,,"Output event list or observation definition file"
  20 +inobs, f, a, "events.fits",,,"Input event list or observation definition file"
  21 +outobs, f, a, "selected_events.fits",,,"Output event list or observation definition file"
22 22 prefix, s, h, "selected_",,,"Prefix for event lists in observation definition file"
  23 +caldb, s, h, "dummy",,, "Calibration database"
  24 +irf, s, h, "cta_dummy_irf",,, "Instrument response function"
23 25 #
24 26 # Selection parameters
25 27 #=====================
... ...
src/ctskymap/ctskymap.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctskymap - Sky mapping tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,6 +35,7 @@
35 35 /* __ Method name definitions ____________________________________________ */
36 36 #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)"
37 37 #define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)"
  38 +#define G_GET_PARAMETERS "ctskymap::get_parameters()"
38 39  
39 40 /* __ Debug definitions __________________________________________________ */
40 41  
... ... @@ -262,12 +263,6 @@ void ctskymap::run(void)
262 263 }
263 264 }
264 265  
265   - // If this is the first valid observation we're working on
266   - // then initialise the sky maps
267   - if (num_obs == 0) {
268   - init_map(obs);
269   - }
270   -
271 266 // Map events into sky map
272 267 map_events(obs);
273 268  
... ... @@ -311,10 +306,10 @@ void ctskymap::save(void)
311 306 }
312 307  
313 308 // Get output filename
314   - m_outfile = (*this)["outfile"].filename();
  309 + m_outmap = (*this)["outmap"].filename();
315 310  
316 311 // Save sky map
317   - m_skymap.save(m_outfile, clobber());
  312 + m_skymap.save(m_outmap, clobber());
318 313  
319 314 // Return
320 315 return;
... ... @@ -332,39 +327,28 @@ void ctskymap::save(void)
332 327 ***************************************************************************/
333 328 void ctskymap::get_parameters(void)
334 329 {
335   - // If there are no observations in container then add a single CTA
336   - // observation using the parameters from the parameter file
  330 + // If there are no observations in container then load them via user
  331 + // parameters
337 332 if (m_obs.size() == 0) {
338 333  
339   - // Get name of CTA events file
340   - m_evfile = (*this)["evfile"].filename();
  334 + // Throw exception if no input observation file is given
  335 + require_inobs(G_GET_PARAMETERS);
341 336  
342   - // Load unbinned CTA observation
343   - GCTAObservation obs;
344   - obs.load(m_evfile);
345   -
346   - // Append CTA observation to container
347   - m_obs.append(obs);
348   -
349   - // Use the xref and yref parameters for binning (otherwise the
350   - // pointing direction(s) is/are used)
351   - m_xref = (*this)["xref"].real();
352   - m_yref = (*this)["yref"].real();
  337 + // Get observation container without response (not needed)
  338 + m_obs = get_observations(false);
353 339  
354 340 } // endif: there was no observation in the container
355 341  
  342 + // Create sky map based on task parameters
  343 + m_skymap = create_map(m_obs);
  344 +
356 345 // Get remaining parameters
357   - m_emin = (*this)["emin"].real();
358   - m_emax = (*this)["emax"].real();
359   - m_proj = (*this)["proj"].string();
360   - m_coordsys = (*this)["coordsys"].string();
361   - m_binsz = (*this)["binsz"].real();
362   - m_nxpix = (*this)["nxpix"].integer();
363   - m_nypix = (*this)["nypix"].integer();
  346 + m_emin = (*this)["emin"].real();
  347 + m_emax = (*this)["emax"].real();
364 348  
365 349 // Read ahead parameters
366 350 if (read_ahead()) {
367   - m_outfile = (*this)["outfile"].filename();
  351 + m_outmap = (*this)["outmap"].filename();
368 352 }
369 353  
370 354 // Return
... ... @@ -373,53 +357,6 @@ void ctskymap::get_parameters(void)
373 357  
374 358  
375 359 /***********************************************************************//**
376   - * @brief Initialise sky map
377   - *
378   - * @param[in] obs CTA observation (no NULL pointer allowed).
379   - *
380   - * This method initialises the sky map.
381   - ***************************************************************************/
382   -void ctskymap::init_map(GCTAObservation* obs)
383   -{
384   - // Clear any existing sky map
385   - m_skymap.clear();
386   -
387   - // Get map centre. If no centre is specified then extract the map centre
388   - // from the pointing direction. This obviously only works for
389   - double xref;
390   - double yref;
391   - if (m_xref != 9999.0 && m_yref != 9999.0) {
392   - xref = m_xref;
393   - yref = m_yref;
394   - }
395   - else {
396   -
397   -
398   - // Get pointer on CTA pointing
399   - const GCTAPointing& pnt = obs->pointing();
400   -
401   - // Set reference point to pointing
402   - if (gammalib::toupper(m_coordsys) == "GAL") {
403   - xref = pnt.dir().l_deg();
404   - yref = pnt.dir().b_deg();
405   - }
406   - else {
407   - xref = pnt.dir().ra_deg();
408   - yref = pnt.dir().dec_deg();
409   - }
410   -
411   - } // endelse: map centre set to pointing
412   -
413   - // Create skymap
414   - m_skymap = GSkymap(m_proj, m_coordsys, xref, yref, -m_binsz, m_binsz,
415   - m_nxpix, m_nypix, 1);
416   -
417   - // Return
418   - return;
419   -}
420   -
421   -
422   -/***********************************************************************//**
423 360 * @brief Map events into a sky map
424 361 *
425 362 * @param[in] obs CTA observation.
... ... @@ -473,8 +410,8 @@ void ctskymap::map_events(GCTAObservation* obs)
473 410 GSkyPixel pixel = m_skymap.dir2pix(dir);
474 411  
475 412 // Skip if pixel is out of range
476   - if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) ||
477   - pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) {
  413 + if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
  414 + pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
478 415 num_outside_map++;
479 416 continue;
480 417 }
... ... @@ -525,19 +462,10 @@ void ctskymap::map_events(GCTAObservation* obs)
525 462 void ctskymap::init_members(void)
526 463 {
527 464 // Initialise members
528   - m_evfile.clear();
529   - m_outfile.clear();
530   - m_proj.clear();
531   - m_coordsys.clear();
532 465 m_obs.clear();
533 466 m_skymap.clear();
534   - m_emin = 0.0;
535   - m_emax = 0.0;
536   - m_xref = 9999.0; // Flags unset (use pointing direction)
537   - m_yref = 9999.0; // Flags unset (use pointing direction)
538   - m_binsz = 0.0;
539   - m_nxpix = 0;
540   - m_nypix = 0;
  467 + m_emin = 0.0;
  468 + m_emax = 0.0;
541 469  
542 470 // Return
543 471 return;
... ... @@ -552,19 +480,10 @@ void ctskymap::init_members(void)
552 480 void ctskymap::copy_members(const ctskymap& app)
553 481 {
554 482 // Copy attributes
555   - m_evfile = app.m_evfile;
556   - m_outfile = app.m_outfile;
557   - m_proj = app.m_proj;
558   - m_coordsys = app.m_coordsys;
559   - m_obs = app.m_obs;
560   - m_skymap = app.m_skymap;
561   - m_emin = app.m_emin;
562   - m_emax = app.m_emax;
563   - m_xref = app.m_xref;
564   - m_yref = app.m_yref;
565   - m_binsz = app.m_binsz;
566   - m_nxpix = app.m_nxpix;
567   - m_nypix = app.m_nypix;
  483 + m_obs = app.m_obs;
  484 + m_skymap = app.m_skymap;
  485 + m_emin = app.m_emin;
  486 + m_emax = app.m_emax;
568 487  
569 488 // Return
570 489 return;
... ...
src/ctskymap/ctskymap.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctskymap - Sky mapping tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -69,21 +69,13 @@ protected:
69 69 void copy_members(const ctskymap& app);
70 70 void free_members(void);
71 71 void get_parameters(void);
72   - void init_map(GCTAObservation* obs);
73 72 void map_events(GCTAObservation* obs);
74 73  
75 74 // User parameters
76   - std::string m_evfile; //!< Input event list
77   - std::string m_outfile; //!< Output counts map
  75 + std::string m_outmap; //!< Output counts map
78 76 double m_emin; //!< Minimum energy (TeV)
79 77 double m_emax; //!< Maximum energy (TeV)
80   - std::string m_proj; //!< WCS projection
81   - std::string m_coordsys; //!< Coordinate system
82   - double m_xref; //!< Longitude reference coordinate
83   - double m_yref; //!< Latitude reference coordinate
84   - double m_binsz; //!< Pixel size
85   - int m_nxpix; //!< Number of pixels in longitude
86   - int m_nypix; //!< Number of pixels in latitude
  78 +
87 79  
88 80 // Protected members
89 81 GObservations m_obs; //!< Observation container
... ...
src/ctskymap/ctskymap.par
  View file @033ed45
... ... @@ -17,8 +17,8 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -evfile, f, a, "events.fits",,,"Event data file name"
21   -outfile, f, a, "skymap.fits",,, "Output file name"
  20 +inobs, f, a, "events.fits",,,"Event data file name"
  21 +outmap, f, a, "skymap.fits",,, "Output file name"
22 22 #
23 23 # Event selection parameters
24 24 #===========================
... ... @@ -27,6 +27,7 @@ emax, r, a, 100.0,,,&quot;Maximum energy in TeV&quot;
27 27 #
28 28 # Spatial binning parameters
29 29 #===========================
  30 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
30 31 nxpix, i, a, 200,,, "Size of the X axis in pixels"
31 32 nypix, i, a, 200,,, "Size of the Y axis in pixels"
32 33 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/cttsmap/cttsmap.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * cttsmap - TS map calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2014-2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -233,8 +233,6 @@ void cttsmap::run(void)
233 233 log.header1("Initialise TS map");
234 234 }
235 235  
236   - // Initialise maps
237   - init_maps();
238 236  
239 237 // Get bins to be computed
240 238 // To split the computation on several jobs,
... ... @@ -369,7 +367,7 @@ void cttsmap::save(void)
369 367 }
370 368  
371 369 // Get output filename
372   - m_outfile = (*this)["outfile"].filename();
  370 + m_outmap = (*this)["outmap"].filename();
373 371  
374 372 // Create fits file instance
375 373 GFits fitsfile = GFits();
... ... @@ -392,7 +390,7 @@ void cttsmap::save(void)
392 390 }
393 391  
394 392 // Save FITS file
395   - fitsfile.saveto(m_outfile,true);
  393 + fitsfile.saveto(m_outmap,true);
396 394  
397 395 // Return
398 396 return;
... ... @@ -411,15 +409,8 @@ void cttsmap::save(void)
411 409 void cttsmap::init_members(void)
412 410 {
413 411 // Initialise members
414   - m_infile.clear();
415   - m_outfile.clear();
416   - m_proj.clear();
417   - m_coordsys.clear();
418   - m_xref = 0.0;
419   - m_yref = 0.0;
420   - m_binsz = 0.0;
421   - m_nxpix = 0;
422   - m_nypix = 0;
  412 + m_srcname.clear();
  413 + m_outmap.clear();
423 414  
424 415 // Initialise protected members
425 416 m_obs.clear();
... ... @@ -445,25 +436,18 @@ void cttsmap::init_members(void)
445 436 void cttsmap::copy_members(const cttsmap& app)
446 437 {
447 438 // Copy attributes
448   - m_infile = app.m_infile;
449   - m_outfile = app.m_outfile;
450   - m_proj = app.m_proj;
451   - m_coordsys = app.m_coordsys;
452   - m_xref = app.m_xref;
453   - m_yref = app.m_yref;
454   - m_binsz = app.m_binsz;
455   - m_nxpix = app.m_nxpix;
456   - m_nypix = app.m_nypix;
  439 + m_srcname = app.m_srcname;
  440 + m_outmap = app.m_outmap;
457 441  
458 442 // Copy protected members
459   - m_binmin = app.m_binmin;
460   - m_binmax = app.m_binmax;
461   - m_logL0 = app.m_logL0;
462   - m_obs = app.m_obs;
463   - m_tsmap = app.m_tsmap;
464   - m_mapnames = app.m_mapnames;
465   - m_maps = app.m_maps;
466   - m_statusmap = app.m_statusmap;
  443 + m_binmin = app.m_binmin;
  444 + m_binmax = app.m_binmax;
  445 + m_logL0 = app.m_logL0;
  446 + m_obs = app.m_obs;
  447 + m_tsmap = app.m_tsmap;
  448 + m_mapnames = app.m_mapnames;
  449 + m_maps = app.m_maps;
  450 + m_statusmap = app.m_statusmap;
467 451  
468 452 // Clone protected members
469 453 m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
... ... @@ -498,42 +482,15 @@ void cttsmap::free_members(void)
498 482 ***************************************************************************/
499 483 void cttsmap::get_parameters(void)
500 484 {
501   - // If there are no observations in container then add a single CTA
502   - // observation using the parameters from the parameter file
  485 + // If there are no observations in container then load them via user
  486 + // parameters
503 487 if (m_obs.size() == 0) {
504 488  
505   - // Allocate CTA observation
506   - GCTAObservation obs;
  489 + // Throw exception if no input observation file is given
  490 + require_inobs(G_GET_PARAMETERS);
507 491  
508   - // Get event file name
509   - std::string filename = (*this)["infile"].filename();
510   -
511   - // Try first to open as FITS file
512   - try {
513   -
514   - // Load data
515   - obs.load(filename);
516   -
517   - // Set response
518   - set_obs_response(&obs);
519   -
520   - // Append observation to container
521   - m_obs.append(obs);
522   -
523   - }
524   -
525   - // ... otherwise try to open as XML file
526   - catch (GException::fits_open_error &e) {
527   -
528   - // Load observations from XML file
529   - m_obs.load(filename);
530   -
531   - // Check if all observations have response information. If
532   - // not, get the calibration database parameters and set
533   - // the response properly
534   - set_response(m_obs);
535   -
536   - } // endcatch: file was an XML file
  492 + // Get observation container
  493 + m_obs = get_observations();
537 494  
538 495 } // endif: there was no observation in the container
539 496  
... ... @@ -542,7 +499,7 @@ void cttsmap::get_parameters(void)
542 499 if (m_obs.models().size() == 0) {
543 500  
544 501 // Get models XML filename
545   - std::string filename = (*this)["srcmdl"].filename();
  502 + std::string filename = (*this)["inmodel"].filename();
546 503  
547 504 // Setup models for optimizing.
548 505 m_obs.models(GModels(filename));
... ... @@ -579,14 +536,11 @@ void cttsmap::get_parameters(void)
579 536 (*m_testsource)["DEC"].fix();
580 537 }
581 538  
582   - // Get map parameters
583   - m_xref = (*this)["xref"].real();
584   - m_yref = (*this)["yref"].real();
585   - m_proj = (*this)["proj"].string();
586   - m_coordsys = (*this)["coordsys"].string();
587   - m_binsz = (*this)["binsz"].real();
588   - m_nxpix = (*this)["nxpix"].integer();
589   - m_nypix = (*this)["nypix"].integer();
  539 + // Create sky map based on task parameters
  540 + GSkymap map = create_map(m_obs);
  541 +
  542 + // Initialise maps from user parameters
  543 + init_maps(map);
590 544  
591 545 // Get optional splitting parameters
592 546 m_binmin = (*this)["binmin"].integer();
... ... @@ -596,7 +550,7 @@ void cttsmap::get_parameters(void)
596 550 // Optionally read ahead parameters so that they get correctly
597 551 // dumped into the log file
598 552 if (read_ahead()) {
599   - m_outfile = (*this)["outfile"].filename();
  553 + m_outmap = (*this)["outmap"].filename();
600 554 }
601 555  
602 556 // Return
... ... @@ -609,23 +563,19 @@ void cttsmap::get_parameters(void)
609 563 *
610 564 * Initialises skymaps that will contain map information.
611 565 ***************************************************************************/
612   -void cttsmap::init_maps(void)
  566 +void cttsmap::init_maps(const GSkymap& map)
613 567 {
614 568 // Initialise map information
615 569 m_tsmap.clear();
616 570  
617 571 // Create skymap
618   - m_tsmap = GSkymap(m_proj, m_coordsys,
619   - m_xref, m_yref, -m_binsz, m_binsz,
620   - m_nxpix, m_nypix, 1);
  572 + m_tsmap = GSkymap(map);
621 573  
622 574 // Initialise map information
623 575 m_statusmap.clear();
624 576  
625 577 // Create status map
626   - m_statusmap = GSkymap(m_proj, m_coordsys,
627   - m_xref, m_yref, -m_binsz, m_binsz,
628   - m_nxpix, m_nypix, 1);
  578 + m_statusmap = GSkymap(map);
629 579  
630 580 // Initialise maps of free parameters
631 581 if (m_testsource != NULL) {
... ... @@ -640,7 +590,7 @@ void cttsmap::init_maps(void)
640 590  
641 591 // Add sky map for free parameters. Note that push_back will
642 592 // create a copy.
643   - m_maps.push_back(m_tsmap);
  593 + m_maps.push_back(map);
644 594  
645 595 // Store parameter name
646 596 m_mapnames.push_back((*m_testsource)[i].name());
... ...
src/cttsmap/cttsmap.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * cttsmap - TS map calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2014-2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTTSMAP_NAME "cttsmap"
37   -#define CTTSMAP_VERSION "00-01-00"
  37 +#define CTTSMAP_VERSION "00-02-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -81,24 +81,16 @@ protected:
81 81 void copy_members(const cttsmap& app);
82 82 void free_members(void);
83 83 void get_parameters(void);
84   - void init_maps(void);
  84 + void init_maps(const GSkymap& map);
85 85  
86 86 // User parameters
87   - std::string m_infile; //!< Input file
88 87 std::string m_srcname; //!< Name of source which is moved around
89   - std::string m_outfile; //!< Output counts map or XML file
90   - std::string m_proj; //!< WCS projection
91   - std::string m_coordsys; //!< Coordinate system
92   - double m_xref; //!< Longitude reference coordinate
93   - double m_yref; //!< Latitude reference coordinate
94   - double m_binsz; //!< Pixel size
95   - int m_nxpix; //!< Number of pixels in longitude
96   - int m_nypix; //!< Number of pixels in latitude
  88 + std::string m_outmap; //!< Output counts map or XML file
97 89  
98 90 // Parameters to control speed and job splitting
99   - int m_binmin; //!< Map bin number from which computation should start
100   - int m_binmax; //!< Map bin number where map computation should end
101   - double m_logL0; //!< Likelihood value of null hypothesis
  91 + int m_binmin; //!< Map bin number from which computation should start
  92 + int m_binmax; //!< Map bin number where map computation should end
  93 + double m_logL0; //!< Likelihood value of null hypothesis
102 94  
103 95 // Protected members
104 96 GObservations m_obs; //!< Observation container
... ...
src/cttsmap/cttsmap.par
  View file @033ed45
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: cttsmap.par #
8   -# Version: 00-01-00 #
  8 +# Version: 00-02-00 #
9 9 # #
10 10 # Author: Michael Mayer #
11 11 # michael.mayer@desy.de #
... ... @@ -17,16 +17,17 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "events.fits",,,"Input event list, cube or observation definition file"
21   -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
22   -srcname, s, a, "TestSource",,,"Test source"
23   -caldb, s, a, "dummy",,, "Calibration database"
24   -irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25   -outfile, f, a, "tsmap.fits",,, "Output counts cube"
  20 +inobs, f, a, "events.fits",,, "Input event list, cube or observation definition file"
  21 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +srcname, s, a, "TestSource",,, "Test source"
  23 +caldb, s, a, "dummy",,, "Calibration database"
  24 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  25 +outmap, f, a, "tsmap.fits",,, "Output counts cube"
26 26  
27 27 #
28 28 # Spatial binning parameters
29 29 #===========================
  30 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
30 31 nxpix, i, a, 200,,, "Size of the X axis in pixels"
31 32 nypix, i, a, 200,,, "Size of the Y axis in pixels"
32 33 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
... ...
src/support/ctool.cpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,7 +47,8 @@
47 47  
48 48  
49 49 /* __ Method name definitions ____________________________________________ */
50   -#define G_GET_EBOUNDS "ctool::get_ebounds()"
  50 +#define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)"
  51 +#define G_CREATE_EBOUNDS "ctool::create_ebounds()"
51 52  
52 53 /* __ Debug definitions __________________________________________________ */
53 54  
... ... @@ -226,6 +227,12 @@ void ctool::init_members(void)
226 227 {
227 228 // Initialise members
228 229 m_read_ahead = false;
  230 + m_use_xml = false;
  231 +
  232 + // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
  233 + // which is defined in GCTALib.hpp. This is somehow a kluge. We need
  234 + // a better mechanism to implement the CTA reference MJD.
  235 + m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
229 236  
230 237 // Set logger properties
231 238 log.date(true);
... ... @@ -244,6 +251,8 @@ void ctool::copy_members(const ctool&amp; app)
244 251 {
245 252 // Copy members
246 253 m_read_ahead = app.m_read_ahead;
  254 + m_use_xml = app.m_use_xml;
  255 + m_cta_ref = app.m_cta_ref;
247 256  
248 257 // Return
249 258 return;
... ... @@ -266,7 +275,102 @@ void ctool::free_members(void)
266 275  
267 276  
268 277 /***********************************************************************//**
269   - * @brief Get the energy boundaries
  278 + * @brief Get observation container
  279 + *
  280 + * @param[in] get_response Indicates whether response information should
  281 + * be loaded (default: true)
  282 + *
  283 + * Get an observation container according to the user parameters. The method
  284 + * supports loading of a individual FITS file or an observation definition
  285 + * file in XML format. If the input filename is empty, parameters are read
  286 + * to build a CTA observation from scratch.
  287 + ***************************************************************************/
  288 +GObservations ctool::get_observations(const bool& get_response)
  289 +{
  290 + // Initialise empty observation container
  291 + GObservations obs;
  292 +
  293 + // Get the filename from the input parameters
  294 + std::string filename = (*this)["inobs"].filename();
  295 +
  296 + // If no observation definition file has been specified then read all
  297 + // parameters that are necessary to create an observation from scratch
  298 + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) {
  299 +
  300 + // Setup a new CTA observation
  301 + GCTAObservation cta_obs = create_cta_obs();
  302 +
  303 + // Get response if required
  304 + if (get_response) {
  305 +
  306 + // Set response
  307 + set_obs_response(&cta_obs);
  308 +
  309 + } // endif: response was required
  310 +
  311 + // Append observation to container
  312 + obs.append(cta_obs);
  313 +
  314 + } // endif: filename was "NONE" or ""
  315 +
  316 + // ... otherwise we have a file name
  317 + else {
  318 +
  319 + // Try first to open as FITS file
  320 + try {
  321 +
  322 + // Allocate empyt CTA observation
  323 + GCTAObservation cta_obs;
  324 +
  325 + // Load data
  326 + cta_obs.load(filename);
  327 +
  328 + // Get response if required
  329 + if (get_response) {
  330 +
  331 + // Set response
  332 + set_obs_response(&cta_obs);
  333 +
  334 + } // endif: response was required
  335 +
  336 + // Append observation to container
  337 + obs.append(cta_obs);
  338 +
  339 + // Signal that no XML file should be used for storage
  340 + m_use_xml = false;
  341 +
  342 + }
  343 +
  344 + // ... otherwise try to open as XML file
  345 + catch (GException::fits_open_error &e) {
  346 +
  347 + // Load observations from XML file
  348 + obs.load(filename);
  349 +
  350 + // Get response if required
  351 + if (get_response) {
  352 +
  353 + // For all observations that have no response, set the response
  354 + // from the task parameters
  355 + set_response(obs);
  356 +
  357 + } // endif: response was required
  358 +
  359 + // Signal that XML file should be used for storage
  360 + m_use_xml = true;
  361 +
  362 + } // endcatch: file was an XML file
  363 +
  364 + }
  365 +
  366 + // Return observation container
  367 + return obs;
  368 +
  369 +}
  370 +
  371 +
  372 +/***********************************************************************//**
  373 + * @brief Create energy boundaries from user parameters
270 374 *
271 375 * @exception GException::invalid_value
272 376 * No valid energy boundary extension found.
... ... @@ -275,8 +379,14 @@ void ctool::free_members(void)
275 379 * supports loading of energy boundary information from the EBOUNDS or
276 380 * ENERGYBINS extension, or setting energy boundaries using a linear
277 381 * or logarithmical spacing.
  382 + *
  383 + * The following parameters are read:
  384 + * ebinalg - Energy binning algorithm
  385 + * emin - Minimum energy (if ebinalg != FILE)
  386 + * emax - Maximum energy (if ebinalg != FILE)
  387 + * enumbins - Number of energy bins (if ebinalg != FILE)
278 388 ***************************************************************************/
279   -GEbounds ctool::get_ebounds(void)
  389 +GEbounds ctool::create_ebounds(void)
280 390 {
281 391 // Allocate energy boundaries
282 392 GEbounds ebounds;
... ... @@ -310,7 +420,7 @@ GEbounds ctool::get_ebounds(void)
310 420 "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
311 421 " is required if the parameter \"ebinalg\""
312 422 " is set to \"FILE\".";
313   - throw GException::invalid_value(G_GET_EBOUNDS, msg);
  423 + throw GException::invalid_value(G_CREATE_EBOUNDS, msg);
314 424 }
315 425  
316 426 } // endif: ebinalg was "FILE"
... ... @@ -343,12 +453,229 @@ GEbounds ctool::get_ebounds(void)
343 453  
344 454  
345 455 /***********************************************************************//**
  456 + * @brief Create a skymap from user parameters
  457 + *
  458 + * @param[in] obs Observation container.
  459 + *
  460 + * Creates an empty skymap from input user parameters. The reference
  461 + * coordinate is extracted from the mean pointing of the CTA observations
  462 + * in the observation container (see get_mean_pointing()).
  463 + *
  464 + * The following parameters are read:
  465 + * usepnt - Use pointing for reference coordinate?
  466 + * xref - Reference in Right Ascension or longitude (if usepnt=no)
  467 + * yref - Reference in Declination or latitude (if usepnt=no)
  468 + * proj - Coordinate projection
  469 + * coordsys - Coordinate system
  470 + * binsz - Bin size (deg)
  471 + * nxpix - Number of pixels in Right Ascension or longitude
  472 + * nypix - Number of pixels in Declination or latitude
  473 + ***************************************************************************/
  474 +GSkymap ctool::create_map(const GObservations& obs)
  475 +{
  476 + // Read task parameters
  477 + double xref = 0.0;
  478 + double yref = 0.0;
  479 + bool usepnt = (*this)["usepnt"].boolean();
  480 + if (!usepnt) {
  481 + xref = (*this)["xref"].real();
  482 + yref = (*this)["yref"].real();
  483 + }
  484 + std::string proj = (*this)["proj"].string();
  485 + std::string coordsys = (*this)["coordsys"].string();
  486 + double binsz = (*this)["binsz"].real();
  487 + int nxpix = (*this)["nxpix"].integer();
  488 + int nypix = (*this)["nypix"].integer();
  489 +
  490 + // If requested, get pointing from observations
  491 + if (usepnt) {
  492 +
  493 + // Get mean pointing
  494 + GSkyDir pnt = get_mean_pointing(obs);
  495 +
  496 + // Use correct coordinate system
  497 + if (gammalib::toupper(coordsys) == "GAL") {
  498 + xref = pnt.l_deg();
  499 + yref = pnt.b_deg();
  500 + }
  501 + else {
  502 + xref = pnt.ra_deg();
  503 + yref = pnt.dec_deg();
  504 + }
  505 +
  506 + } // endif: got mean pointing as reference
  507 +
  508 + // Initialise sky map
  509 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz,
  510 + nxpix, nypix, 1);
  511 +
  512 + // Return sky map
  513 + return map;
  514 +}
  515 +
  516 +
  517 +/***********************************************************************//**
  518 + * @brief Create a CTA event cube from user parameters
  519 + *
  520 + * @param[in] obs Observation container.
  521 + *
  522 + * Creates an empty CTA event cube from input user parameters. The reference
  523 + * coordinate is extracted from the mean pointing of the CTA observations
  524 + * in the observation container (see get_mean_pointing()). This method will
  525 + * be used if a ctool has the input parameter usepnt=true. The method
  526 + * appends a dummy GTI to the event cube.
  527 + *
  528 + * The following parameters are read:
  529 + * proj - Coordinate projection
  530 + * coordsys - Coordinate system
  531 + * binsz - Bin size (deg)
  532 + * nxpix - Number of pixels in Right Ascension or longitude
  533 + * nypix - Number of pixels in Declination or latitude
  534 + * ebinalg - Energy binning algorithm
  535 + * emin - Minimum energy (if ebinalg != FILE)
  536 + * emax - Maximum energy (if ebinalg != FILE)
  537 + * enumbins - Number of energy bins (if ebinalg != FILE)
  538 + ***************************************************************************/
  539 +GCTAEventCube ctool::create_cube(const GObservations& obs)
  540 +{
  541 + // Get skymap
  542 + GSkymap map = create_map(obs);
  543 +
  544 + // Set energy boundaries
  545 + GEbounds ebounds = create_ebounds();
  546 +
  547 + // Set dummy GTI that is needed for event cube creation
  548 + GGti gti;
  549 + gti.append(GTime(0.0), GTime(0.1234));
  550 +
  551 + // Extend skymap to the requested number of maps
  552 + map.nmaps(ebounds.size());
  553 +
  554 + // Initialise event cube from all parameters
  555 + GCTAEventCube cube = GCTAEventCube(map, ebounds, gti);
  556 +
  557 + // Return event cube
  558 + return cube;
  559 +}
  560 +
  561 +
  562 +/***********************************************************************//**
  563 + * @brief Create a CTA observation from user parameters
  564 + *
  565 + * Creates an empty CTA observation from user parameters. An empty event list
  566 + * including RoI, GTI and energy boundary information is attached to the
  567 + * observation. The method also sets the pointing direction using the ra and
  568 + * dec parameter, the ROI based on ra, dec and rad, a single GTI based on
  569 + * tmin and tmax, and a single energy boundary based on emin and emax. The
  570 + * method furthermore sets the ontime, livetime and deadtime correction
  571 + * factor.
  572 + *
  573 + * The following parameters are read:
  574 + * ra - Right Ascension of pointing and RoI centre (deg)
  575 + * dec - Declination of pointing and RoI centre (deg)
  576 + * rad - Radius of RoI (deg)
  577 + * deadc - Deadtime correction factor
  578 + * tmin - Start time
  579 + * tmax - Stop time
  580 + * emin - Minimum energy (TeV)
  581 + * emax - Maximum energy (TeV)
  582 + ***************************************************************************/
  583 +GCTAObservation ctool::create_cta_obs(void)
  584 +{
  585 + // Get CTA observation parameters
  586 + double ra = (*this)["ra"].real();
  587 + double dec = (*this)["dec"].real();
  588 + double rad = (*this)["rad"].real();
  589 + double deadc = (*this)["deadc"].real();
  590 + double t_min = (*this)["tmin"].real();
  591 + double t_max = (*this)["tmax"].real();
  592 + double e_min = (*this)["emin"].real();
  593 + double e_max = (*this)["emax"].real();
  594 +
  595 + // Allocate CTA observation and empty event list
  596 + GCTAObservation obs;
  597 + GCTAEventList list;
  598 +
  599 + // Set pointing direction
  600 + GCTAPointing pnt;
  601 + GSkyDir skydir;
  602 + skydir.radec_deg(ra, dec);
  603 + pnt.dir(skydir);
  604 +
  605 + // Set ROI
  606 + GCTAInstDir instdir(skydir);
  607 + GCTARoi roi(instdir, rad);
  608 +
  609 + // Set GTI
  610 + GGti gti(m_cta_ref);
  611 + GTime tstart;
  612 + GTime tstop;
  613 + tstart.set(t_min, m_cta_ref);
  614 + tstop.set(t_max, m_cta_ref);
  615 + gti.append(tstart, tstop);
  616 +
  617 + // Set energy boundaries
  618 + GEbounds ebounds;
  619 + GEnergy emin;
  620 + GEnergy emax;
  621 + emin.TeV(e_min);
  622 + emax.TeV(e_max);
  623 + ebounds.append(emin, emax);
  624 +
  625 + // Set CTA event list attributes
  626 + list.roi(roi);
  627 + list.gti(gti);
  628 + list.ebounds(ebounds);
  629 +
  630 + // Attach event list to CTA observation
  631 + obs.events(list);
  632 +
  633 + // Set observation ontime, livetime and deadtime correction factor
  634 + obs.pointing(pnt);
  635 + obs.ontime(gti.ontime());
  636 + obs.livetime(gti.ontime()*deadc);
  637 + obs.deadc(deadc);
  638 +
  639 + // Return CTA observation
  640 + return obs;
  641 +
  642 +}
  643 +
  644 +
  645 +/***********************************************************************//**
  646 + * @brief Throws exception in inobs parameter is not valid
  647 + *
  648 + * @param[in] method Method name.
  649 + ***************************************************************************/
  650 +void ctool::require_inobs(const std::string& method)
  651 +{
  652 + // Throw exception if no infile is given
  653 + if ((*this)["inobs"].filename() == "NONE" ||
  654 + (*this)["inobs"].filename() == "") {
  655 + std::string msg = "A valid file needs to be specified for the "
  656 + "\"inobs\" parameter, yet \""+
  657 + (*this)["inobs"].filename()+"\" was given."
  658 + " Specify a valid observation definition or "
  659 + "FITS file to proceed.";
  660 + throw GException::invalid_value(method, msg);
  661 + }
  662 +
  663 + // Return
  664 + return;
  665 +}
  666 +
  667 +
  668 +/***********************************************************************//**
346 669 * @brief Set response for all CTA observations in container
347 670 *
348 671 * @param[in,out] obs Observation container
349 672 *
350 673 * Set the response for a CTA observations in the container that so far have
351   - * no response using the "database" and "irf" task parameters.
  674 + * no response using the "caldb" and "irf" task parameters.
  675 + *
  676 + * The following parameters are read:
  677 + * caldb - Calibration database
  678 + * irf - Instrument response function
352 679 ***************************************************************************/
353 680 void ctool::set_response(GObservations& obs)
354 681 {
... ... @@ -380,8 +707,12 @@ void ctool::set_response(GObservations&amp; obs)
380 707 *
381 708 * @param[in,out] obs CTA observation
382 709 *
383   - * Set the response for a CTA observation using the "database" and "irf"
  710 + * Set the response for a CTA observation using the "caldb" and "irf"
384 711 * task parameters.
  712 + *
  713 + * The following parameters are read:
  714 + * caldb - Calibration database
  715 + * irf - Instrument response function
385 716 ***************************************************************************/
386 717 void ctool::set_obs_response(GCTAObservation* obs)
387 718 {
... ... @@ -409,6 +740,67 @@ void ctool::set_obs_response(GCTAObservation* obs)
409 740  
410 741  
411 742 /***********************************************************************//**
  743 + * @brief Derives mean pointing from CTA observations
  744 + *
  745 + * @param[in] obs Observation container.
  746 + * @return Pointing direction.
  747 + *
  748 + * Computes the mean pointing direction from all CTA observations in the
  749 + * observation container.
  750 + *
  751 + * @todo This method does not work properly for wrap arounds in Right
  752 + * Ascension. It certainly will also lead to bad results near the
  753 + * celestial pole.
  754 + ***************************************************************************/
  755 +GSkyDir ctool::get_mean_pointing(const GObservations& obs)
  756 +{
  757 + // Initialise values
  758 + double xref = 0.0;
  759 + double yref = 0.0;
  760 + int n_pnt = 0;
  761 +
  762 + // Loop over all observations to compute mean pointings
  763 + for (int i = 0; i < obs.size(); ++i) {
  764 +
  765 + // If we have a CTA observation then extract the pointing direction
  766 + const GCTAObservation* cta_obs = dynamic_cast<const GCTAObservation*>(obs[i]);
  767 + if (cta_obs != NULL) {
  768 +
  769 + // Get pointing
  770 + const GCTAPointing& pnt = cta_obs->pointing();
  771 +
  772 + // Add to coordinates
  773 + xref += pnt.dir().ra_deg();
  774 + yref += pnt.dir().dec_deg();
  775 + n_pnt += 1;
  776 +
  777 + } // endif: cta observation was valid
  778 +
  779 + } // endfor: loop over all observations
  780 +
  781 + // Signal if no pointing is found
  782 + if (n_pnt < 1) {
  783 + std::string msg = "No valid CTA observation has been found in "
  784 + "observation list, hence no pointing information "
  785 + "could be extracted. Use the \"usepnt=no\" "
  786 + "option and specify the pointing explicitly.";
  787 + throw GException::invalid_value(G_GET_MEAN_POINTING, msg);
  788 + }
  789 +
  790 + // Calculate mean coordinates
  791 + double ra = xref / double(n_pnt);
  792 + double dec = yref / double(n_pnt);
  793 +
  794 + // Initialise sky dir
  795 + GSkyDir dir = GSkyDir();
  796 + dir.radec_deg(ra,dec);
  797 +
  798 + // Return sky direction
  799 + return dir;
  800 +}
  801 +
  802 +
  803 +/***********************************************************************//**
412 804 * @brief Get current resident set size (physical memory use) in Bytes
413 805 *
414 806 * @return Physical memory use in Bytes.
... ...
src/support/ctool.hpp
  View file @033ed45
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -60,19 +60,37 @@ public:
60 60 // Public methods
61 61 virtual void execute(void);
62 62  
  63 +#ifdef SWIG
  64 +public:
  65 +#else
63 66 protected:
  67 +#endif
64 68 // Protected methods
65   - void init_members(void);
66   - void copy_members(const ctool& app);
67   - void free_members(void);
68   - const bool& read_ahead(void) const;
69   - GEbounds get_ebounds(void);
70   - void set_response(GObservations& obs);
71   - void set_obs_response(GCTAObservation* obs);
72   - size_t get_current_rss(void);
  69 + void init_members(void);
  70 + void copy_members(const ctool& app);
  71 + void free_members(void);
  72 + const bool& read_ahead(void) const;
  73 + GObservations get_observations(const bool& get_response = true);
  74 +
  75 + // Protected methods that create objects from user parameters
  76 + GEbounds create_ebounds(void);
  77 + GSkymap create_map(const GObservations& obs);
  78 + GCTAEventCube create_cube(const GObservations& obs);
  79 + GCTAObservation create_cta_obs(void);
  80 +
  81 + // Protected methods that check user parameters
  82 + void require_inobs(const std::string& method);
  83 +
  84 + // Protected support methods
  85 + void set_response(GObservations& obs);
  86 + void set_obs_response(GCTAObservation* obs);
  87 + GSkyDir get_mean_pointing(const GObservations& obs);
  88 + size_t get_current_rss(void);
73 89  
74 90 // Protected members
75   - bool m_read_ahead; //!< Read ahead parameters
  91 + bool m_read_ahead; //!< Read ahead parameters
  92 + bool m_use_xml; //!< Use XML file instead of FITS file for observations
  93 + GTimeReference m_cta_ref; //!< CTA time reference
76 94 };
77 95  
78 96  
... ...
test/test_cscripts.sh.in
  View file @033ed45
... ... @@ -52,7 +52,7 @@ export PFILES=pfiles
52 52 $ECHO @ECHO_N@ "Test cspull: "
53 53 chmod +x @abs_top_srcdir@/scripts/cspull.py
54 54 @abs_top_srcdir@/scripts/cspull.py \
55   - srcmdl="data/crab.xml" \
  55 + inmodel="data/crab.xml" \
56 56 outfile="pull.dat" \
57 57 ntrials="3" \
58 58 caldb="irf" \
... ... @@ -138,9 +138,9 @@ $ECHO &quot; ok&quot;
138 138 $ECHO @ECHO_N@ "Test csresmap: "
139 139 chmod +x @abs_top_srcdir@/scripts/csresmap.py
140 140 @abs_top_srcdir@/scripts/csresmap.py \
141   - infile="data/crab_events.fits" \
142   - outfile="resmap.fits" \
143   - srcmdl="data/crab.xml" \
  141 + inobs="data/crab_events.fits" \
  142 + outmap="resmap.fits" \
  143 + inmodel="data/crab.xml" \
144 144 caldb="irf" \
145 145 irf="cta_dummy_irf" \
146 146 nxpix=50 \
... ...
test/test_ctbin.py
  View file @033ed45
... ... @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 64 """
65 65 # Set-up ctbin
66 66 bin = ctools.ctbin()
67   - bin["evfile"].filename(self.events_name)
68   - bin["outfile"].filename("cntmap.fits")
  67 + bin["inobs"].filename(self.events_name)
  68 + bin["outcube"].filename("cntmap.fits")
69 69 bin["ebinalg"].string("LOG")
70 70 bin["emin"].real(0.1)
71 71 bin["emax"].real(100.0)
... ...
test/test_ctbkgcube.py
  View file @033ed45
... ... @@ -40,6 +40,8 @@ class Test(gammalib.GPythonTestSuite):
40 40 # Set members
41 41 self.events_name = "data/crab_events.fits"
42 42 self.bkg_model = "data/crab.xml"
  43 + self.caldb = "irf"
  44 + self.irf = "cta_dummy_irf"
43 45  
44 46 # Return
45 47 return
... ... @@ -65,11 +67,13 @@ class Test(gammalib.GPythonTestSuite):
65 67 """
66 68 # Set-up ctbkgcube
67 69 bkgcube = ctools.ctbkgcube()
68   - bkgcube["infile"].filename(self.events_name)
  70 + bkgcube["inobs"].filename(self.events_name)
69 71 bkgcube["inmodel"].filename(self.bkg_model)
70   - bkgcube["cntmap"].filename("NONE")
71   - bkgcube["outfile"].filename("bkgcube.fits")
  72 + bkgcube["incube"].filename("NONE")
  73 + bkgcube["outcube"].filename("bkgcube.fits")
72 74 bkgcube["outmodel"].filename("bkgcube.xml")
  75 + bkgcube["caldb"].string(self.caldb)
  76 + bkgcube["irf"].string(self.irf)
73 77 bkgcube["ebinalg"].string("LOG")
74 78 bkgcube["emin"].real(0.1)
75 79 bkgcube["emax"].real(100.0)
... ...
test/test_ctbutterfly.py
  View file @033ed45
... ... @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite):
67 67 """
68 68 # Set-up ctbutterfly
69 69 butterfly = ctools.ctbutterfly()
70   - butterfly["infile"].filename(self.events_name)
71   - butterfly["srcmdl"].filename(self.model_name)
  70 + butterfly["inobs"].filename(self.events_name)
  71 + butterfly["inmodel"].filename(self.model_name)
72 72 butterfly["srcname"].string("Crab")
73 73 butterfly["caldb"].string(self.caldb)
74 74 butterfly["irf"].string(self.irf)
... ...
test/test_ctcubemask.py
  View file @033ed45
... ... @@ -65,9 +65,9 @@ class Test(gammalib.GPythonTestSuite):
65 65 """
66 66 # Set-up ctcubemask
67 67 mask = ctools.ctcubemask()
68   - mask["infile"].filename(self.cntmap_name)
  68 + mask["inobs"].filename(self.cntmap_name)
69 69 mask["regfile"].filename(self.regfile)
70   - mask["outfile"].filename("filtered_cntmap.fits")
  70 + mask["outcube"].filename("filtered_cntmap.fits")
71 71 mask["ra"].real(83.63)
72 72 mask["dec"].real(22.01)
73 73 mask["rad"].real(2.0)
... ...
test/test_ctexpcube.py
  View file @033ed45
... ... @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite):
66 66 """
67 67 # Set-up ctexpcube
68 68 expcube = ctools.ctexpcube()
69   - expcube["infile"].filename(self.events_name)
70   - expcube["cntmap"].filename("NONE")
71   - expcube["outfile"].filename("expcube.fits")
  69 + expcube["inobs"].filename(self.events_name)
  70 + expcube["incube"].filename("NONE")
  71 + expcube["outcube"].filename("expcube.fits")
72 72 expcube["caldb"].string(self.caldb)
73 73 expcube["irf"].string(self.irf)
74 74 expcube["ebinalg"].string("LOG")
... ...
test/test_ctlike.py
  View file @033ed45
... ... @@ -67,11 +67,11 @@ class Test(gammalib.GPythonTestSuite):
67 67 """
68 68 # Set-up ctlike
69 69 like = ctools.ctlike()
70   - like["infile"].filename(self.events_name)
71   - like["srcmdl"].filename(self.model_name)
  70 + like["inobs"].filename(self.events_name)
  71 + like["inmodel"].filename(self.model_name)
72 72 like["caldb"].string(self.caldb)
73 73 like["irf"].string(self.irf)
74   - like["outmdl"].filename("result.xml")
  74 + like["outmodel"].filename("result.xml")
75 75  
76 76 # Run tool
77 77 self.test_try("Run ctlike")
... ...
test/test_ctmodel.py
  View file @033ed45
... ... @@ -66,12 +66,13 @@ class Test(gammalib.GPythonTestSuite):
66 66 """
67 67 # Set-up ctmodel from scratch
68 68 model = ctools.ctmodel()
69   - model["infile"].filename("NONE")
70   - model["outfile"].filename("modmap.fits")
71   - model["srcmdl"].filename(self.model_name)
72   - model["obsfile"].filename("NONE")
  69 + model["incube"].filename("NONE")
  70 + model["outcube"].filename("modmap.fits")
  71 + model["inmodel"].filename(self.model_name)
  72 + model["inobs"].filename("NONE")
73 73 model["caldb"].string(self.caldb)
74 74 model["irf"].string(self.irf)
  75 + model["rad"].real(5.0)
75 76 model["ra"].real(83.63)
76 77 model["dec"].real(22.01)
77 78 model["tmin"].real(0.0)
... ...
test/test_ctobssim.py
  View file @033ed45
... ... @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite):
67 67 """
68 68 # Set-up ctobssim
69 69 sim = ctools.ctobssim()
70   - sim["infile"].filename(self.model_name)
71   - sim["outfile"].filename("events.fits")
  70 + sim["inmodel"].filename(self.model_name)
  71 + sim["outevents"].filename("events.fits")
72 72 sim["caldb"].string(self.caldb)
73 73 sim["irf"].string(self.irf)
74 74 sim["ra"].real(83.63)
... ... @@ -117,8 +117,8 @@ class Test(gammalib.GPythonTestSuite):
117 117  
118 118 # Set-up ctobssim
119 119 sim = ctools.ctobssim(obs)
120   - sim["outfile"].filename("sim_events.xml")
121   -
  120 + sim["outevents"].filename("sim_events.xml")
  121 + sim["inmodel"].filename(self.model_name)
122 122 # Run tool
123 123 self.test_try("Run ctobssim")
124 124 try:
... ...
test/test_ctools.sh
  View file @033ed45
... ... @@ -67,10 +67,11 @@ export PFILES=pfiles
67 67 # Test ctobssim
68 68 # =============
69 69 $ECHO -n "Test ctobssim: "
  70 +
70 71 #
71 72 # Run 1
72   -$ctobssim infile="data/crab.xml" \
73   - outfile="events.fits" \
  73 +$ctobssim inmodel="data/crab.xml" \
  74 + outevents="events.fits" \
74 75 caldb="irf" \
75 76 irf="cta_dummy_irf" \
76 77 ra=83.63 \
... ... @@ -90,8 +91,8 @@ else
90 91 fi
91 92 #
92 93 # Run 2
93   -$ctobssim infile="data/model_background_cube.xml" \
94   - outfile="events2.fits" \
  94 +$ctobssim inmodel="data/model_background_cube.xml" \
  95 + outevents="events2.fits" \
95 96 caldb="irf" \
96 97 irf="cta_dummy_irf" \
97 98 ra=84.17263 \
... ... @@ -116,8 +117,8 @@ $ECHO &quot; ok&quot;
116 117 # Test ctskymap
117 118 # =============
118 119 $ECHO -n "Test ctskymap: "
119   -$ctskymap evfile="events.fits" \
120   - outfile="skymap.fits" \
  120 +$ctskymap inobs="events.fits" \
  121 + outmap="skymap.fits" \
121 122 emin=0.1 \
122 123 emax=100.0 \
123 124 nxpix=200 \
... ... @@ -144,8 +145,8 @@ $ECHO &quot; ok&quot;
144 145 $ECHO -n "Test ctbin: "
145 146 #
146 147 # Run 1
147   -$ctbin evfile="events.fits" \
148   - outfile="cntmap.fits" \
  148 +$ctbin inobs="events.fits" \
  149 + outcube="cntmap.fits" \
149 150 emin=0.1 \
150 151 emax=100.0 \
151 152 enumbins=20 \
... ... @@ -167,8 +168,8 @@ else
167 168 fi
168 169 #
169 170 # Run 2
170   -$ctbin evfile="events2.fits" \
171   - outfile="cntmap2.fits" \
  171 +$ctbin inobs="events2.fits" \
  172 + outcube="cntmap2.fits" \
172 173 emin=0.3 \
173 174 emax=8.0 \
174 175 enumbins=10 \
... ... @@ -191,18 +192,19 @@ fi
191 192 $ECHO " ok"
192 193  
193 194  
  195 +
194 196 #
195 197 # Test ctmodel
196 198 # ============
197 199 $ECHO -n "Test ctmodel: "
198 200 #
199 201 # Run 1
200   -$ctmodel infile="cntmap.fits" \
201   - obsfile="cntmap.fits" \
202   - outfile="modmap1.fits" \
  202 +$ctmodel inobs="cntmap.fits" \
  203 + incube="cntmap.fits" \
  204 + outcube="modmap1.fits" \
203 205 caldb="irf" \
204 206 irf="cta_dummy_irf" \
205   - srcmdl="data/crab.xml"
  207 + inmodel="data/crab.xml"
206 208 $ECHO -n "."
207 209 if [ -s "modmap1.fits" ]
208 210 then
... ... @@ -213,14 +215,15 @@ else
213 215 fi
214 216 #
215 217 # Run 3
216   -$ctmodel infile="NONE" \
217   - obsfile="NONE" \
218   - outfile="modmap3.fits" \
  218 +$ctmodel inobs="NONE" \
  219 + incube="NONE" \
  220 + outcube="modmap3.fits" \
219 221 caldb="irf" \
220 222 irf="cta_dummy_irf" \
221   - srcmdl="data/crab.xml" \
  223 + inmodel="data/crab.xml" \
222 224 ra=83.63 \
223 225 dec=22.01 \
  226 + rad=10.0 \
224 227 tmin=0.0 \
225 228 tmax=1800.0 \
226 229 emin=0.1 \
... ... @@ -243,12 +246,12 @@ else
243 246 fi
244 247 #
245 248 # Run 2
246   -$ctmodel infile="cntmap2.fits" \
247   - obsfile="cntmap2.fits" \
248   - outfile="modmap2.fits" \
  249 +$ctmodel inobs="cntmap2.fits" \
  250 + incube="cntmap2.fits" \
  251 + outcube="modmap2.fits" \
249 252 caldb="irf" \
250 253 irf="cta_dummy_irf" \
251   - srcmdl="data/model_background_cube.xml"
  254 + inmodel="data/model_background_cube.xml"
252 255 $ECHO -n "."
253 256 if [ -s "modmap2.fits" ]
254 257 then
... ... @@ -264,8 +267,8 @@ $ECHO &quot; ok&quot;
264 267 # Test ctselect
265 268 # =============
266 269 $ECHO -n "Test ctselect: "
267   -$ctselect infile="events.fits" \
268   - outfile="selected_events.fits" \
  270 +$ctselect inobs="events.fits" \
  271 + outobs="selected_events.fits" \
269 272 ra=83.63 \
270 273 dec=22.01 \
271 274 rad=10.0 \
... ... @@ -290,9 +293,9 @@ $ECHO &quot; ok&quot;
290 293 $ECHO -n "Test ctlike: "
291 294 #
292 295 # Run 1
293   -$ctlike infile="cntmap.fits" \
294   - srcmdl="data/crab.xml" \
295   - outmdl="results_binned.xml" \
  296 +$ctlike inobs="cntmap.fits" \
  297 + inmodel="data/crab.xml" \
  298 + outmodel="results_binned.xml" \
296 299 caldb="irf" \
297 300 irf="cta_dummy_irf"
298 301 $ECHO -n "."
... ... @@ -305,9 +308,9 @@ else
305 308 fi
306 309 #
307 310 # Run 2
308   -$ctlike infile="selected_events.fits" \
309   - srcmdl="data/crab.xml" \
310   - outmdl="results_unbinned.xml" \
  311 +$ctlike inobs="selected_events.fits" \
  312 + inmodel="data/crab.xml" \
  313 + outmodel="results_unbinned.xml" \
311 314 caldb="irf" \
312 315 irf="cta_dummy_irf"
313 316 $ECHO -n "."
... ... @@ -320,9 +323,9 @@ else
320 323 fi
321 324 #
322 325 # Run 3
323   -$ctlike infile="events2.fits" \
324   - srcmdl="data/model_background_cube.xml" \
325   - outmdl="results_unbinned_background_cube.xml" \
  326 +$ctlike inobs="events2.fits" \
  327 + inmodel="data/model_background_cube.xml" \
  328 + outmodel="results_unbinned_background_cube.xml" \
326 329 caldb="irf" \
327 330 irf="cta_dummy_irf"
328 331 $ECHO -n "."
... ... @@ -335,9 +338,9 @@ else
335 338 fi
336 339 #
337 340 # Run 4
338   -$ctlike infile="cntmap2.fits" \
339   - srcmdl="data/model_background_cube.xml" \
340   - outmdl="results_binned_background_cube.xml" \
  341 +$ctlike inobs="cntmap2.fits" \
  342 + inmodel="data/model_background_cube.xml" \
  343 + outmodel="results_binned_background_cube.xml" \
341 344 caldb="irf" \
342 345 irf="cta_dummy_irf"
343 346 $ECHO -n "."
... ... @@ -355,12 +358,12 @@ $ECHO &quot; ok&quot;
355 358 # Test cttsmap
356 359 # ============
357 360 $ECHO -n "Test cttsmap: "
358   -$cttsmap infile="selected_events.fits" \
359   - srcmdl="data/crab.xml" \
  361 +$cttsmap inobs="selected_events.fits" \
  362 + inmodel="data/crab.xml" \
360 363 srcname="Crab" \
361 364 caldb="irf" \
362 365 irf="cta_dummy_irf" \
363   - outfile="tsmap.fits" \
  366 + outmap="tsmap.fits" \
364 367 nxpix=5 \
365 368 nypix=5 \
366 369 binsz=0.05 \
... ... @@ -385,9 +388,9 @@ $ECHO &quot; ok&quot;
385 388 $ECHO -n "Test ctexpcube: "
386 389 #
387 390 # Run 1
388   -$ctexpcube infile="data/crab_events.fits" \
389   - cntmap="NONE" \
390   - outfile="expcube1.fits" \
  391 +$ctexpcube inobs="data/crab_events.fits" \
  392 + incube="NONE" \
  393 + outcube="expcube1.fits" \
391 394 caldb="irf" \
392 395 irf="cta_dummy_irf" \
393 396 emin=0.1 \
... ... @@ -410,9 +413,9 @@ else
410 413 fi
411 414 #
412 415 # Run 2
413   -$ctexpcube infile="data/crab_events.fits" \
414   - cntmap="data/crab_cntmap.fits" \
415   - outfile="expcube2.fits" \
  416 +$ctexpcube inobs="data/crab_events.fits" \
  417 + incube="data/crab_cntmap.fits" \
  418 + outcube="expcube2.fits" \
416 419 caldb="irf" \
417 420 irf="cta_dummy_irf"
418 421 $ECHO -n "."
... ... @@ -432,9 +435,9 @@ $ECHO &quot; ok&quot;
432 435 $ECHO -n "Test ctpsfcube: "
433 436 #
434 437 # Run 1
435   -$ctpsfcube infile="data/crab_events.fits" \
436   - cntmap="NONE" \
437   - outfile="psfcube1.fits" \
  438 +$ctpsfcube inobs="data/crab_events.fits" \
  439 + incube="NONE" \
  440 + outcube="psfcube1.fits" \
438 441 caldb="irf" \
439 442 irf="cta_dummy_irf" \
440 443 emin=0.1 \
... ... @@ -459,9 +462,9 @@ else
459 462 fi
460 463 #
461 464 # Run 2
462   -$ctpsfcube infile="data/crab_events.fits" \
463   - cntmap="data/crab_cntmap.fits" \
464   - outfile="psfcube2.fits" \
  465 +$ctpsfcube inobs="data/crab_events.fits" \
  466 + incube="data/crab_cntmap.fits" \
  467 + outcube="psfcube2.fits" \
465 468 caldb="irf" \
466 469 irf="cta_dummy_irf" \
467 470 amax=0.3 \
... ... @@ -483,11 +486,13 @@ $ECHO &quot; ok&quot;
483 486 $ECHO -n "Test ctbkgcube: "
484 487 #
485 488 # Run 1
486   -$ctbkgcube infile="data/crab_events.fits" \
  489 +$ctbkgcube inobs="data/crab_events.fits" \
487 490 inmodel="data/crab.xml" \
488   - cntmap="NONE" \
489   - outfile="bkgcube1.fits" \
  491 + incube="NONE" \
  492 + outcube="bkgcube1.fits" \
490 493 outmodel="bkgcube1.xml" \
  494 + caldb="irf" \
  495 + irf="cta_dummy_irf" \
491 496 emin=0.1 \
492 497 emax=100.0 \
493 498 enumbins=20 \
... ... @@ -515,10 +520,12 @@ else
515 520 fi
516 521 #
517 522 # Run 2
518   -$ctbkgcube infile="data/crab_events.fits" \
  523 +$ctbkgcube inobs="data/crab_events.fits" \
519 524 inmodel="data/crab.xml" \
520   - cntmap="data/crab_cntmap.fits" \
521   - outfile="bkgcube2.fits" \
  525 + incube="data/crab_cntmap.fits" \
  526 + caldb="irf" \
  527 + irf="cta_dummy_irf" \
  528 + outcube="bkgcube2.fits" \
522 529 outmodel="bkgcube2.xml"
523 530 $ECHO -n "."
524 531 if [ -s "bkgcube2.fits" ]
... ... @@ -542,9 +549,9 @@ $ECHO &quot; ok&quot;
542 549 # Test ctcubemask
543 550 # ===============
544 551 $ECHO -n "Test ctcubemask: "
545   -$ctcubemask infile="data/crab_cntmap.fits" \
  552 +$ctcubemask inobs="data/crab_cntmap.fits" \
546 553 regfile="data/exclusion.reg" \
547   - outfile="filtered_cube.fits" \
  554 + outcube="filtered_cube.fits" \
548 555 ra=83.63 \
549 556 dec=22.01 \
550 557 rad=2.0 \
... ... @@ -565,8 +572,8 @@ $ECHO &quot; ok&quot;
565 572 # Test ctbutterfly
566 573 # ===============
567 574 $ECHO -n "Test ctbutterfly: "
568   -$ctbutterfly infile="data/crab_events.fits.gz" \
569   - srcmdl="data/crab.xml" \
  575 +$ctbutterfly inobs="data/crab_events.fits.gz" \
  576 + inmodel="data/crab.xml" \
570 577 srcname="Crab" \
571 578 outfile="butterfly.txt" \
572 579 caldb="irf" \
... ...
test/test_ctpsfcube.py
  View file @033ed45
... ... @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite):
66 66 """
67 67 # Set-up ctpsfcube
68 68 psfcube = ctools.ctpsfcube()
69   - psfcube["infile"].filename(self.events_name)
70   - psfcube["cntmap"].filename("NONE")
71   - psfcube["outfile"].filename("psfcube.fits")
  69 + psfcube["inobs"].filename(self.events_name)
  70 + psfcube["incube"].filename("NONE")
  71 + psfcube["outcube"].filename("psfcube.fits")
72 72 psfcube["caldb"].string(self.caldb)
73 73 psfcube["irf"].string(self.irf)
74 74 psfcube["ebinalg"].string("LOG")
... ...
test/test_ctselect.py
  View file @033ed45
... ... @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 64 """
65 65 # Set-up ctselect
66 66 select = ctools.ctselect()
67   - select["infile"].filename(self.events_name)
68   - select["outfile"].filename("selected_events.fits")
  67 + select["inobs"].filename(self.events_name)
  68 + select["outobs"].filename("selected_events.fits")
69 69 select["ra"].real(83.63)
70 70 select["dec"].real(22.01)
71 71 select["rad"].real(3.0)
... ...
test/test_ctskymap.py
  View file @033ed45
... ... @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 64 """
65 65 # Set-up ctskymap
66 66 skymap = ctools.ctskymap()
67   - skymap["evfile"].filename(self.events_name)
68   - skymap["outfile"].filename("skymap.fits")
  67 + skymap["inobs"].filename(self.events_name)
  68 + skymap["outmap"].filename("skymap.fits")
69 69 skymap["emin"].real(0.1)
70 70 skymap["emax"].real(100.0)
71 71 skymap["nxpix"].integer(200)
... ...
test/test_cttsmap.py
  View file @033ed45
... ... @@ -67,12 +67,12 @@ class Test(gammalib.GPythonTestSuite):
67 67 """
68 68 # Set-up cttsmap
69 69 tsmap = ctools.cttsmap()
70   - tsmap["infile"].filename(self.events_name)
71   - tsmap["srcmdl"].filename(self.model_name)
  70 + tsmap["inobs"].filename(self.events_name)
  71 + tsmap["inmodel"].filename(self.model_name)
72 72 tsmap["srcname"].string("Crab")
73 73 tsmap["caldb"].string(self.caldb)
74 74 tsmap["irf"].string(self.irf)
75   - tsmap["outfile"].filename("tsmap.fits")
  75 + tsmap["outmap"].filename("tsmap.fits")
76 76 tsmap["nxpix"].integer(5)
77 77 tsmap["nypix"].integer(5)
78 78 tsmap["binsz"].real(0.02)
... ...
test/test_pipelines.py
  View file @033ed45
... ... @@ -78,8 +78,8 @@ class Test(gammalib.GPythonTestSuite):
78 78  
79 79 # Simulate events
80 80 sim = ctools.ctobssim()
81   - sim["infile"].filename(model_name)
82   - sim["outfile"].filename(events_name)
  81 + sim["inmodel"].filename(model_name)
  82 + sim["outevents"].filename(events_name)
83 83 sim["caldb"].string(caldb)
84 84 sim["irf"].string(irf)
85 85 sim["ra"].real(ra)
... ... @@ -98,8 +98,8 @@ class Test(gammalib.GPythonTestSuite):
98 98  
99 99 # Select events
100 100 select = ctools.ctselect()
101   - select["infile"].filename(events_name)
102   - select["outfile"].filename(selected_events_name)
  101 + select["inobs"].filename(events_name)
  102 + select["outobs"].filename(selected_events_name)
103 103 select["ra"].real(ra)
104 104 select["dec"].real(dec)
105 105 select["rad"].real(rad_select)
... ... @@ -116,9 +116,9 @@ class Test(gammalib.GPythonTestSuite):
116 116  
117 117 # Perform maximum likelihood fitting
118 118 like = ctools.ctlike()
119   - like["infile"].filename(selected_events_name)
120   - like["srcmdl"].filename(model_name)
121   - like["outmdl"].filename(result_name)
  119 + like["inobs"].filename(selected_events_name)
  120 + like["inmodel"].filename(model_name)
  121 + like["outmodel"].filename(result_name)
122 122 like["caldb"].string(caldb)
123 123 like["irf"].string(irf)
124 124 self.test_try("Execute ctlike")
... ... @@ -149,7 +149,7 @@ class Test(gammalib.GPythonTestSuite):
149 149  
150 150 # Simulate events
151 151 sim = ctools.ctobssim()
152   - sim["infile"].filename(model_name)
  152 + sim["inmodel"].filename(model_name)
153 153 sim["caldb"].string(caldb)
154 154 sim["irf"].string(irf)
155 155 sim["ra"].real(ra)
... ...
test/test_python.py
  View file @033ed45
... ... @@ -102,6 +102,7 @@ if __name__ == &#39;__main__&#39;:
102 102  
103 103 # Run test suite
104 104 success = suites.run()
  105 +
105 106  
106 107 # Save test results
107 108 suites.save("reports/ctools.xml")
... ...