Commit be3e4fa07e18a76816e99234bd10930217919f5e

Authored by Michael Mayer
1 parent 70701752

initial commit to extend ctool-base class and unify/simplify the get_parameters(…

…)-function of individual ctools
pyext/ctools/obsutils.py
@@ -42,6 +42,7 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, @@ -42,6 +42,7 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0,
42 binsz - Pixel size for binned simulation (deg/pixel) 42 binsz - Pixel size for binned simulation (deg/pixel)
43 npix - Number of pixels in X and Y for binned simulation 43 npix - Number of pixels in X and Y for binned simulation
44 """ 44 """
  45 +
45 # Allocate ctobssim application and set parameters 46 # Allocate ctobssim application and set parameters
46 sim = ctools.ctobssim(obs) 47 sim = ctools.ctobssim(obs)
47 sim["seed"].integer(seed) 48 sim["seed"].integer(seed)
scripts/cspull.py
@@ -19,10 +19,12 @@ @@ -19,10 +19,12 @@
19 # 19 #
20 # ========================================================================== 20 # ==========================================================================
21 import gammalib 21 import gammalib
  22 +import ctools
  23 +#import obsutils
22 from ctools import obsutils 24 from ctools import obsutils
23 import sys 25 import sys
24 import csv 26 import csv
25 -import math 27 +#import math
26 28
27 29
28 # ============ # 30 # ============ #
@@ -41,12 +43,12 @@ class cspull(gammalib.GApplication): @@ -41,12 +43,12 @@ class cspull(gammalib.GApplication):
41 """ 43 """
42 # Set name 44 # Set name
43 self.name = "cspull" 45 self.name = "cspull"
44 - self.version = "0.3.0" 46 + self.version = "0.4.0"
45 47
46 # Initialise some members 48 # Initialise some members
47 self.obs = None 49 self.obs = None
48 self.model = None 50 self.model = None
49 - self.m_srcmdl = None 51 + self.m_inmodel = None
50 52
51 # Make sure that parfile exists 53 # Make sure that parfile exists
52 file = self.parfile() 54 file = self.parfile()
@@ -93,7 +95,7 @@ class cspull(gammalib.GApplication): @@ -93,7 +95,7 @@ class cspull(gammalib.GApplication):
93 95
94 # Create default parfile 96 # Create default parfile
95 pars = gammalib.GApplicationPars() 97 pars = gammalib.GApplicationPars()
96 - pars.append(gammalib.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model")) 98 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
97 pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name")) 99 pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
98 pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) 100 pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
99 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database")) 101 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
@@ -121,9 +123,18 @@ class cspull(gammalib.GApplication): @@ -121,9 +123,18 @@ class cspull(gammalib.GApplication):
121 """ 123 """
122 Get parameters from parfile and setup the observation. 124 Get parameters from parfile and setup the observation.
123 """ 125 """
124 - # Get parameters 126 +
  127 +# sim = ctools.ctobssim()
  128 +# sim.run()
  129 +# self.obs = sim.obs()
  130 +# sim = ctools.ctobssim()
  131 +# sim.run()
  132 +# self.obs = sim.obs().copy()
  133 +# print self.obs
  134 +#
  135 +# # Get parameters
125 if self.model == None: 136 if self.model == None:
126 - self.m_srcmdl = self["srcmdl"].filename() 137 + self.m_inmodel = self["inmodel"].filename()
127 self.m_outfile = self["outfile"].filename() 138 self.m_outfile = self["outfile"].filename()
128 self.m_ntrials = self["ntrials"].integer() 139 self.m_ntrials = self["ntrials"].integer()
129 self.m_caldb = self["caldb"].string() 140 self.m_caldb = self["caldb"].string()
@@ -144,15 +155,15 @@ class cspull(gammalib.GApplication): @@ -144,15 +155,15 @@ class cspull(gammalib.GApplication):
144 155
145 # Set some fixed parameters 156 # Set some fixed parameters
146 self.m_log = False # Logging in client tools 157 self.m_log = False # Logging in client tools
147 - self.m_debug = False # Debugging in client tools 158 + self.m_debug = True # Debugging in client tools
148 159
149 # Setup observations 160 # Setup observations
150 self.obs = self.set_obs() 161 self.obs = self.set_obs()
151 - 162 +
152 # Load source model 163 # Load source model
153 - if self.m_srcmdl != None:  
154 - self.model = gammalib.GModels(self.m_srcmdl)  
155 - 164 + if self.m_inmodel != None:
  165 + self.model = gammalib.GModels(self.m_inmodel)
  166 +
156 # Append source model to observation 167 # Append source model to observation
157 self.obs.models(self.model) 168 self.obs.models(self.model)
158 169
@@ -260,18 +271,23 @@ class cspull(gammalib.GApplication): @@ -260,18 +271,23 @@ class cspull(gammalib.GApplication):
260 if self.logExplicit(): 271 if self.logExplicit():
261 self.log.header2("Trial "+str(seed+1)) 272 self.log.header2("Trial "+str(seed+1))
262 273
  274 +# sim = ctools.ctobssim(self.obs)
  275 +# sim["seed"].integer(seed)
  276 +# sim.run()
  277 + print self.obs
263 # Simulate events 278 # Simulate events
264 - obs = obsutils.sim(self.obs, \  
265 - nbins=self.m_enumbins, \  
266 - seed=seed, \  
267 - binsz=self.m_binsz, \  
268 - npix=self.m_npix, \  
269 - edisp=self.m_edisp, \  
270 - log=self.m_log, debug=self.m_debug)  
271 - 279 + sim_obs = obsutils.sim(self.obs, \
  280 + nbins=self.m_enumbins, \
  281 + seed=seed, \
  282 + binsz=self.m_binsz, \
  283 + npix=self.m_npix, \
  284 + edisp=self.m_edisp, \
  285 + log=self.m_log, debug=self.m_debug)
  286 + print "SIM"
  287 + print sim_obs
272 # Determine number of events in simulation 288 # Determine number of events in simulation
273 nevents = 0.0 289 nevents = 0.0
274 - for run in obs: 290 + for run in sim_obs:
275 nevents += run.events().number() 291 nevents += run.events().number()
276 292
277 # Write simulation results 293 # Write simulation results
@@ -282,8 +298,10 @@ class cspull(gammalib.GApplication): @@ -282,8 +298,10 @@ class cspull(gammalib.GApplication):
282 self.log("\n") 298 self.log("\n")
283 299
284 # Fit model 300 # Fit model
285 - like = obsutils.fit(obs, edisp=self.m_edisp, \  
286 - log=self.m_log, debug=self.m_debug) 301 + like = ctools.ctlike(sim_obs)
  302 + like.run()
  303 +# like = obsutils.fit(obs, edisp=self.m_edisp, \
  304 +# log=self.m_log, debug=self.m_debug)
287 305
288 # Store results 306 # Store results
289 logL = like.opt().value() 307 logL = like.opt().value()
@@ -320,7 +338,7 @@ class cspull(gammalib.GApplication): @@ -320,7 +338,7 @@ class cspull(gammalib.GApplication):
320 338
321 # Compute pull 339 # Compute pull
322 fitted_value = par.value() 340 fitted_value = par.value()
323 - real_value = self.model[i][k].value() 341 + real_value = self.obs.models()[i][k].value()
324 error = par.error() 342 error = par.error()
325 if error != 0.0: 343 if error != 0.0:
326 pull = (fitted_value - real_value) / error 344 pull = (fitted_value - real_value) / error
scripts/csresmap.py
@@ -21,10 +21,7 @@ @@ -21,10 +21,7 @@
21 import ctools 21 import ctools
22 import gammalib 22 import gammalib
23 import sys 23 import sys
24 -import tempfile  
25 import os 24 import os
26 -import math  
27 -  
28 25
29 # ============== # 26 # ============== #
30 # cstsdist class # 27 # cstsdist class #
@@ -42,7 +39,7 @@ class csresmap(gammalib.GApplication): @@ -42,7 +39,7 @@ class csresmap(gammalib.GApplication):
42 """ 39 """
43 # Set name 40 # Set name
44 self.name = "csresmap" 41 self.name = "csresmap"
45 - self.version = "0.1.0" 42 + self.version = "0.2.0"
46 43
47 # Initialise some members 44 # Initialise some members
48 self.obs = None 45 self.obs = None
@@ -103,9 +100,9 @@ class csresmap(gammalib.GApplication): @@ -103,9 +100,9 @@ class csresmap(gammalib.GApplication):
103 100
104 # Create default parfile 101 # Create default parfile
105 pars = gammalib.GApplicationPars() 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 pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name")) 106 pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
110 pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name")) 107 pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
111 pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) 108 pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
@@ -137,7 +134,7 @@ class csresmap(gammalib.GApplication): @@ -137,7 +134,7 @@ class csresmap(gammalib.GApplication):
137 # Set observation if not done before 134 # Set observation if not done before
138 if self.obs.size() == 0: 135 if self.obs.size() == 0:
139 136
140 - obsfile = self["infile"].filename() 137 + obsfile = self["inobs"].filename()
141 try: 138 try:
142 139
143 self.obs = gammalib.GObservations(obsfile) 140 self.obs = gammalib.GObservations(obsfile)
@@ -163,8 +160,8 @@ class csresmap(gammalib.GApplication): @@ -163,8 +160,8 @@ class csresmap(gammalib.GApplication):
163 # Note that the models should be optimised for 160 # Note that the models should be optimised for
164 # meaningful residual map 161 # meaningful residual map
165 if self.obs.models().size() == 0: 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 self.m_modfile = self["modfile"].filename() 165 self.m_modfile = self["modfile"].filename()
169 self.m_binfile = self["binfile"].filename() 166 self.m_binfile = self["binfile"].filename()
170 self.m_xref = self["xref"].real() 167 self.m_xref = self["xref"].real()
scripts/cstsdist.py
@@ -346,6 +346,7 @@ class cstsdist(GApplication): @@ -346,6 +346,7 @@ class cstsdist(GApplication):
346 346
347 # Set source name 347 # Set source name
348 source.name("Test") 348 source.name("Test")
  349 + source.tscalc(True)
349 350
350 # Return source 351 # Return source
351 return source 352 return source
@@ -368,7 +369,7 @@ class cstsdist(GApplication): @@ -368,7 +369,7 @@ class cstsdist(GApplication):
368 binsz=self.m_binsz, \ 369 binsz=self.m_binsz, \
369 npix=self.m_npix, \ 370 npix=self.m_npix, \
370 log=self.m_log, debug=self.m_debug) 371 log=self.m_log, debug=self.m_debug)
371 - 372 +
372 # Determine number of events in simulation 373 # Determine number of events in simulation
373 nevents = 0.0 374 nevents = 0.0
374 for run in sim: 375 for run in sim:
@@ -381,6 +382,11 @@ class cstsdist(GApplication): @@ -381,6 +382,11 @@ class cstsdist(GApplication):
381 self.log(nevents) 382 self.log(nevents)
382 self.log("\n") 383 self.log("\n")
383 384
  385 +# sim.models(self.full_model)
  386 +
  387 +# like = obsutils.fit(sim, log = self.m_log, debug = self.m_debug)
  388 +# ts = like.obs().models()["Test"].ts()
  389 +
384 # Fit background only 390 # Fit background only
385 sim.models(self.bkg_model) 391 sim.models(self.bkg_model)
386 like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) 392 like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
@@ -407,6 +413,7 @@ class cstsdist(GApplication): @@ -407,6 +413,7 @@ class cstsdist(GApplication):
407 # Fit background and test source 413 # Fit background and test source
408 sim.models(self.full_model) 414 sim.models(self.full_model)
409 like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) 415 like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  416 + print like_all.obs().models()
410 result_all = like_all.obs().models() 417 result_all = like_all.obs().models()
411 LogL_all = like_all.opt().value() 418 LogL_all = like_all.opt().value()
412 npred_all = like_all.obs().npred() 419 npred_all = like_all.obs().npred()
src/ctbin/ctbin.cpp
@@ -33,8 +33,8 @@ @@ -33,8 +33,8 @@
33 #include "GTools.hpp" 33 #include "GTools.hpp"
34 34
35 /* __ Method name definitions ____________________________________________ */ 35 /* __ Method name definitions ____________________________________________ */
36 -#define G_INIT_CUBE "ctbin::init_cube()"  
37 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)" 36 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)"
  37 +#define G_GET_PARAMETERS "ctbin::get_parameters()"
38 38
39 /* __ Debug definitions __________________________________________________ */ 39 /* __ Debug definitions __________________________________________________ */
40 40
@@ -238,9 +238,6 @@ void ctbin::run(void) @@ -238,9 +238,6 @@ void ctbin::run(void)
238 } 238 }
239 } 239 }
240 240
241 - // Initialise counts cube  
242 - init_cube();  
243 -  
244 // Loop over all observations in the container 241 // Loop over all observations in the container
245 for (int i = 0; i < m_obs.size(); ++i) { 242 for (int i = 0; i < m_obs.size(); ++i) {
246 243
@@ -309,14 +306,14 @@ void ctbin::save(void) @@ -309,14 +306,14 @@ void ctbin::save(void)
309 } 306 }
310 307
311 // Get output filename 308 // Get output filename
312 - m_outfile = (*this)["outfile"].filename(); 309 + m_outcube = (*this)["outcube"].filename();
313 310
314 // Get CTA observation from observation container 311 // Get CTA observation from observation container
315 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]); 312 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
316 313
317 // Save only if observation is valid 314 // Save only if observation is valid
318 if (obs != NULL) { 315 if (obs != NULL) {
319 - obs->save(m_outfile, clobber()); 316 + obs->save(m_outcube, clobber());
320 } 317 }
321 318
322 // Return 319 // Return
@@ -336,16 +333,8 @@ void ctbin::save(void) @@ -336,16 +333,8 @@ void ctbin::save(void)
336 void ctbin::init_members(void) 333 void ctbin::init_members(void)
337 { 334 {
338 // Initialise members 335 // Initialise members
339 - m_evfile.clear();  
340 - m_outfile.clear();  
341 - m_proj.clear();  
342 - m_coordsys.clear(); 336 + m_outcube.clear();
343 m_usepnt = false; 337 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;  
349 338
350 // Initialise protected members 339 // Initialise protected members
351 m_obs.clear(); 340 m_obs.clear();
@@ -368,16 +357,8 @@ void ctbin::init_members(void) @@ -368,16 +357,8 @@ void ctbin::init_members(void)
368 void ctbin::copy_members(const ctbin& app) 357 void ctbin::copy_members(const ctbin& app)
369 { 358 {
370 // Copy attributes 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; 360 + m_outcube = app.m_outcube;
375 m_usepnt = app.m_usepnt; 361 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;  
381 362
382 // Copy protected members 363 // Copy protected members
383 m_obs = app.m_obs; 364 m_obs = app.m_obs;
@@ -405,6 +386,9 @@ void ctbin::free_members(void) @@ -405,6 +386,9 @@ void ctbin::free_members(void)
405 /***********************************************************************//** 386 /***********************************************************************//**
406 * @brief Get application parameters 387 * @brief Get application parameters
407 * 388 *
  389 + * @exception GException::invalid_value
  390 + * Parameter "inobs" is required for ctbin.
  391 + *
408 * Get all task parameters from parameter file or (if required) by querying 392 * Get all task parameters from parameter file or (if required) by querying
409 * the user. Most parameters are only required if no observation exists so 393 * the user. Most parameters are only required if no observation exists so
410 * far in the observation container. In this case, a single CTA observation 394 * far in the observation container. In this case, a single CTA observation
@@ -413,131 +397,60 @@ void ctbin::free_members(void) @@ -413,131 +397,60 @@ void ctbin::free_members(void)
413 ***************************************************************************/ 397 ***************************************************************************/
414 void ctbin::get_parameters(void) 398 void ctbin::get_parameters(void)
415 { 399 {
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 400 + // If there are no observations in container then load them via user parameters
421 if (m_obs.size() == 0) { 401 if (m_obs.size() == 0) {
422 402
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 {  
431 -  
432 - // Load event list in CTA observation  
433 - obs.load(m_evfile);  
434 -  
435 - // Append CTA observation to container  
436 - m_obs.append(obs); 403 + // Throw exception if no infile is given, since this tool needs an observation
  404 + // including events
  405 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
437 406
  407 + std::string msg = "Parameter \"inobs\" is required to be given in ctbin."
  408 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  409 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
438 } 410 }
439 411
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);  
445 -  
446 - } 412 + // Build observation container without response (not needed)
  413 + m_obs = get_observations(false);
447 414
448 } // endif: there was no observation in the container 415 } // endif: there was no observation in the container
449 416
450 // Get parameters 417 // Get parameters
451 m_usepnt = (*this)["usepnt"].boolean(); 418 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();  
464 419
465 - // Optionally read ahead parameters so that they get correctly  
466 - // dumped into the log file  
467 - if (read_ahead()) {  
468 - m_outfile = (*this)["outfile"].filename();  
469 - } 420 + // Initialise event cube
  421 + GCTAEventCube cube;
470 422
471 - // Return  
472 - return;  
473 -} 423 + // Check if pointing should be used
  424 + if (m_usepnt) {
474 425
  426 + // buld cube from user parameters and pointing information
  427 + cube = get_cube(m_obs);
475 428
476 -/***********************************************************************//**  
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(); 429 + } // endif: pointing used as map center
492 430
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) { 431 + else {
498 432
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 - } 433 + // Build cube from user parameters
  434 + cube = get_cube();
518 435
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 - } 436 + } // endelse: m_usepnt was false
528 437
529 - } // endif: used pointing 438 + // Get the skymap
  439 + m_cube = cube.map();
530 440
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()); 441 + // Get energy boundaries
  442 + m_ebounds = cube.ebounds();
  443 +
  444 + // Optionally read ahead parameters so that they get correctly
  445 + // dumped into the log file
  446 + if (read_ahead()) {
  447 + m_outcube = (*this)["outcube"].filename();
  448 + }
535 449
536 // Return 450 // Return
537 return; 451 return;
538 } 452 }
539 453
540 -  
541 /***********************************************************************//** 454 /***********************************************************************//**
542 * @brief Fill events into counts cube 455 * @brief Fill events into counts cube
543 * 456 *
@@ -590,8 +503,8 @@ void ctbin::fill_cube(GCTAObservation* obs) @@ -590,8 +503,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
590 } 503 }
591 504
592 // Skip if pixel is out of range 505 // 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)) { 506 + if (pixel.x() < -0.5 || pixel.x() > (m_cube.nx()-0.5) ||
  507 + pixel.y() < -0.5 || pixel.y() > (m_cube.ny()-0.5)) {
595 num_outside_map++; 508 num_outside_map++;
596 continue; 509 continue;
597 } 510 }
src/ctbin/ctbin.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTBIN_NAME "ctbin" 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,21 +87,12 @@ protected:
87 void copy_members(const ctbin& app); 87 void copy_members(const ctbin& app);
88 void free_members(void); 88 void free_members(void);
89 void get_parameters(void); 89 void get_parameters(void);
90 - void init_cube(void);  
91 void fill_cube(GCTAObservation* obs); 90 void fill_cube(GCTAObservation* obs);
92 void obs_cube(void); 91 void obs_cube(void);
93 92
94 // User parameters 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 bool m_usepnt; //!< Use pointing instead of xref/yref parameters 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 // Protected members 97 // Protected members
107 GObservations m_obs; //!< Observation container 98 GObservations m_obs; //!< Observation container
src/ctbin/ctbin.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctbin.par # 7 # File: ctbin.par #
8 -# Version: 00-05-00 # 8 +# Version: 00-06-00 #
9 # # 9 # #
10 # Author: Jurgen Knodlseder # 10 # Author: Jurgen Knodlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -17,8 +17,8 @@ @@ -17,8 +17,8 @@
17 # 17 #
18 # General parameters 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 # Energy binning parameters 23 # Energy binning parameters
24 #========================== 24 #==========================
src/ctbkgcube/ctbkgcube.cpp
@@ -34,8 +34,6 @@ @@ -34,8 +34,6 @@
34 34
35 /* __ Method name definitions ____________________________________________ */ 35 /* __ Method name definitions ____________________________________________ */
36 #define G_GET_PARAMETERS "ctbkgcube::get_parameters()" 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 #define G_FILL_CUBE "ctbkgcube::fill_cube(GCTAObservation*)" 37 #define G_FILL_CUBE "ctbkgcube::fill_cube(GCTAObservation*)"
40 38
41 /* __ Debug definitions __________________________________________________ */ 39 /* __ Debug definitions __________________________________________________ */
@@ -407,7 +405,7 @@ void ctbkgcube::save(void) @@ -407,7 +405,7 @@ void ctbkgcube::save(void)
407 } 405 }
408 406
409 // Get output filenames 407 // Get output filenames
410 - std::string outfile = (*this)["outfile"].filename(); 408 + std::string outfile = (*this)["outcube"].filename();
411 std::string outmodel = (*this)["outmodel"].filename(); 409 std::string outmodel = (*this)["outmodel"].filename();
412 410
413 // Create empty FITS file 411 // Create empty FITS file
@@ -460,6 +458,8 @@ void ctbkgcube::save(void) @@ -460,6 +458,8 @@ void ctbkgcube::save(void)
460 void ctbkgcube::init_members(void) 458 void ctbkgcube::init_members(void)
461 { 459 {
462 // Initialise members 460 // Initialise members
  461 + m_outcube.clear();
  462 + m_outmodel.clear();
463 m_obs.clear(); 463 m_obs.clear();
464 m_bkgcube.clear(); 464 m_bkgcube.clear();
465 m_bkgmdl.clear(); 465 m_bkgmdl.clear();
@@ -481,6 +481,8 @@ void ctbkgcube::init_members(void) @@ -481,6 +481,8 @@ void ctbkgcube::init_members(void)
481 void ctbkgcube::copy_members(const ctbkgcube& app) 481 void ctbkgcube::copy_members(const ctbkgcube& app)
482 { 482 {
483 // Copy members 483 // Copy members
  484 + m_outmodel = app.m_outmodel;
  485 + m_outcube = app.m_outcube;
484 m_obs = app.m_obs; 486 m_obs = app.m_obs;
485 m_bkgcube = app.m_bkgcube; 487 m_bkgcube = app.m_bkgcube;
486 m_bkgmdl = app.m_bkgmdl; 488 m_bkgmdl = app.m_bkgmdl;
@@ -507,6 +509,8 @@ void ctbkgcube::free_members(void) @@ -507,6 +509,8 @@ void ctbkgcube::free_members(void)
507 /***********************************************************************//** 509 /***********************************************************************//**
508 * @brief Get application parameters 510 * @brief Get application parameters
509 * 511 *
  512 + * @exception GException::invalid_value
  513 + * Parameter "inobs" is required for ctbkgcube.
510 * @exception GException::invalid_value 514 * @exception GException::invalid_value
511 * No background model definition XML file specified. 515 * No background model definition XML file specified.
512 * 516 *
@@ -515,11 +519,23 @@ void ctbkgcube::free_members(void) @@ -515,11 +519,23 @@ void ctbkgcube::free_members(void)
515 ***************************************************************************/ 519 ***************************************************************************/
516 void ctbkgcube::get_parameters(void) 520 void ctbkgcube::get_parameters(void)
517 { 521 {
518 - // If we do not have any observations in the container then load them  
519 - // using the "infile" parameter 522 + // If there are no observations in container then load them via user parameters
520 if (m_obs.size() == 0) { 523 if (m_obs.size() == 0) {
521 - get_obs();  
522 - } 524 +
  525 + // Throw exception if no infile is given, since this tool needs an observation
  526 + // including events
  527 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
  528 +
  529 + std::string msg = "Parameter \"inobs\" is required to be given in ctbkgcube."
  530 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  531 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  532 + }
  533 +
  534 + // Build observation container
  535 + m_obs = get_observations();
  536 +
  537 + } // endif: there was no observation in the container
  538 +
523 539
524 // If there are no models associated with the observations then load now 540 // If there are no models associated with the observations then load now
525 // the model definition from the XML file 541 // the model definition from the XML file
@@ -536,176 +552,54 @@ void ctbkgcube::get_parameters(void) @@ -536,176 +552,54 @@ void ctbkgcube::get_parameters(void)
536 m_obs.models(models); 552 m_obs.models(models);
537 } 553 }
538 554
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); 555 + // Get the incube filename
  556 + std::string incube = (*this)["incube"].filename();
566 557
567 - } 558 + // Check for filename validity
  559 + if ((gammalib::toupper(incube) == "NONE") ||
  560 + (gammalib::strip_whitespace(incube) == "")) {
568 561
569 - // ... otherwise setup the background cube from the counts map  
570 - else {  
571 -  
572 - // Set background cube from counts map  
573 - set_from_cntmap(cntmap);  
574 -  
575 - } 562 + // Check for usepnt flag
  563 + bool usepnt = (*this)["usepnt"].boolean();
576 564
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 - } 565 + // Check if pointing should be used
  566 + if (usepnt) {
610 567
611 - // Try first to open as FITS file  
612 - try { 568 + // build cube from user parameters and pointing information
  569 + m_bkgcube = get_cube(m_obs);
613 570
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 -} 571 + } // endif: pointing used as map center
636 572
  573 + else {
637 574
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); 575 + // Build cube from user parameters
  576 + m_bkgcube = get_cube();
653 577
654 - // Set background cube from counts map  
655 - const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events()); 578 + } // endelse: m_usepnt was false
656 579
657 - // Continue only if cube is valid  
658 - if (cube != NULL) { 580 + } // endif: filename was not valid
659 581
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 - } 582 + // ... otherwise setup the background cube from the counts map
  583 + else {
  584 +
  585 + // Set background cube from counts map
  586 + m_bkgcube = set_from_cntmap(incube);
  587 +
  588 + } // endelse: cube was loaded from file
694 589
695 - } // endif: observation contained an events cube 590 + // Get energy definition
  591 + m_ebounds = m_bkgcube.ebounds();
696 592
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); 593 + // Read output filenames (if needed)
  594 + if (read_ahead()) {
  595 + m_outcube = (*this)["outcube"].filename();
  596 + m_outmodel = (*this)["outmodel"].filename();
702 } 597 }
703 598
704 // Return 599 // Return
705 return; 600 return;
706 } 601 }
707 602
708 -  
709 /***********************************************************************//** 603 /***********************************************************************//**
710 * @brief Generate background cube 604 * @brief Generate background cube
711 * 605 *
src/ctbkgcube/ctbkgcube.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTBKGCUBE_NAME "ctbkgcube" 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,12 @@ protected: @@ -68,11 +68,12 @@ protected:
68 void copy_members(const ctbkgcube& app); 68 void copy_members(const ctbkgcube& app);
69 void free_members(void); 69 void free_members(void);
70 void get_parameters(void); 70 void get_parameters(void);
71 - void get_obs(void);  
72 - void set_from_cntmap(const std::string& filename);  
73 void fill_cube(GCTAObservation* obs); 71 void fill_cube(GCTAObservation* obs);
74 72
75 // Protected members 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 GObservations m_obs; //!< Observation container
77 GCTAEventCube m_bkgcube; //!< Background cube 78 GCTAEventCube m_bkgcube; //!< Background cube
78 GModels m_bkgmdl; //!< CTA background models 79 GModels m_bkgmdl; //!< CTA background models
src/ctbkgcube/ctbkgcube.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctbkgcube.par # 7 # File: ctbkgcube.par #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Chia-Chun Lu # 10 # Author: Chia-Chun Lu #
11 # c.lu@mpi-hd.mpg.de # 11 # c.lu@mpi-hd.mpg.de #
@@ -15,13 +15,15 @@ @@ -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 inmodel, f, a, "NONE",,, "Input (background) model XML file" 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 outmodel, f, a, "NONE",,,"Output (background) model XML file" 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 # Energy binning parameters 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,6 +35,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
33 # 35 #
34 # Spatial binning parameters 36 # Spatial binning parameters
35 #=========================== 37 #===========================
  38 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
36 nxpix, i, a, 200,,, "Size of the X axis in pixels" 39 nxpix, i, a, 200,,, "Size of the X axis in pixels"
37 nypix, i, a, 200,,, "Size of the Y axis in pixels" 40 nypix, i, a, 200,,, "Size of the Y axis in pixels"
38 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 41 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/ctbutterfly/ctbutterfly.cpp
@@ -399,7 +399,7 @@ void ctbutterfly::save(void) @@ -399,7 +399,7 @@ void ctbutterfly::save(void)
399 void ctbutterfly::init_members(void) 399 void ctbutterfly::init_members(void)
400 { 400 {
401 // Initialise members 401 // Initialise members
402 - m_infile.clear(); 402 + //m_inobs.clear();
403 m_srcname.clear(); 403 m_srcname.clear();
404 m_outfile.clear(); 404 m_outfile.clear();
405 m_ebounds.clear(); 405 m_ebounds.clear();
@@ -424,7 +424,7 @@ void ctbutterfly::init_members(void) @@ -424,7 +424,7 @@ void ctbutterfly::init_members(void)
424 void ctbutterfly::copy_members(const ctbutterfly& app) 424 void ctbutterfly::copy_members(const ctbutterfly& app)
425 { 425 {
426 // Copy attributes 426 // Copy attributes
427 - m_infile = app.m_infile; 427 + // m_inobs = app.m_inobs;
428 m_srcname = app.m_srcname; 428 m_srcname = app.m_srcname;
429 m_outfile = app.m_outfile; 429 m_outfile = app.m_outfile;
430 m_ebounds = app.m_ebounds; 430 m_ebounds = app.m_ebounds;
@@ -455,6 +455,9 @@ void ctbutterfly::free_members(void) @@ -455,6 +455,9 @@ void ctbutterfly::free_members(void)
455 * @brief Get application parameters 455 * @brief Get application parameters
456 * 456 *
457 * @exception GException::invalid_value 457 * @exception GException::invalid_value
  458 + * Parameter "inobs" is required for ctbutterfly.
  459 +
  460 + * @exception GException::invalid_value
458 * Test source not found or no RA/DEC parameters found for test 461 * Test source not found or no RA/DEC parameters found for test
459 * source. 462 * source.
460 * 463 *
@@ -466,51 +469,30 @@ void ctbutterfly::free_members(void) @@ -466,51 +469,30 @@ void ctbutterfly::free_members(void)
466 ***************************************************************************/ 469 ***************************************************************************/
467 void ctbutterfly::get_parameters(void) 470 void ctbutterfly::get_parameters(void)
468 { 471 {
469 - // If there are no observations in container then add a single CTA  
470 - // observation using the parameters from the parameter file 472 + // If there are no observations in container then load them via user parameters
471 if (m_obs.size() == 0) { 473 if (m_obs.size() == 0) {
472 474
473 - // Allocate CTA observation  
474 - GCTAObservation obs;  
475 -  
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); 475 + // Throw exception if no infile is given, since this tool needs an observation
  476 + // including events
  477 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
490 478
  479 + std::string msg = "Parameter \"inobs\" is required to be given in ctbutterfly."
  480 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  481 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
491 } 482 }
492 483
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 484 + // Build observation container
  485 + m_obs = get_observations();
505 486
506 } // endif: there was no observation in the container 487 } // endif: there was no observation in the container
507 488
  489 +
508 // If there is are no models associated with the observations then 490 // If there is are no models associated with the observations then
509 // load now the model definition 491 // load now the model definition
510 if (m_obs.models().size() == 0) { 492 if (m_obs.models().size() == 0) {
511 493
512 // Get models XML filename 494 // Get models XML filename
513 - std::string filename = (*this)["srcmdl"].filename(); 495 + std::string filename = (*this)["inmodel"].filename();
514 496
515 // Setup models for optimizing. 497 // Setup models for optimizing.
516 m_obs.models(GModels(filename)); 498 m_obs.models(GModels(filename));
src/ctbutterfly/ctbutterfly.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTBUTTERFLY_NAME "ctbutterfly" 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,7 +80,6 @@ protected: @@ -80,7 +80,6 @@ protected:
80 void get_parameters(void); 80 void get_parameters(void);
81 81
82 // User parameters 82 // User parameters
83 - std::string m_infile; //!< Input file  
84 std::string m_srcname; //!< Name of source to compute butterfly 83 std::string m_srcname; //!< Name of source to compute butterfly
85 std::string m_outfile; //!< Output ascii file 84 std::string m_outfile; //!< Output ascii file
86 GEbounds m_ebounds; //!< Energy binning definition 85 GEbounds m_ebounds; //!< Energy binning definition
src/ctbutterfly/ctbutterfly.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctbutterfly # 7 # File: ctbutterfly #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Michael Mayer # 10 # Author: Michael Mayer #
11 # michael.mayer@physik.hu-berlin.de # 11 # michael.mayer@physik.hu-berlin.de #
@@ -17,8 +17,8 @@ @@ -17,8 +17,8 @@
17 # 17 #
18 # General parameters 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 srcname, s, a, "Crab",,,"Source of interest" 22 srcname, s, a, "Crab",,,"Source of interest"
23 caldb, s, a, "dummy",,, "Calibration database" 23 caldb, s, a, "dummy",,, "Calibration database"
24 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
src/ctcubemask/ctcubemask.cpp
@@ -36,6 +36,7 @@ @@ -36,6 +36,7 @@
36 /* __ Method name definitions ____________________________________________ */ 36 /* __ Method name definitions ____________________________________________ */
37 #define G_RUN "ctcubemask::run()" 37 #define G_RUN "ctcubemask::run()"
38 #define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)" 38 #define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)"
  39 +#define G_GET_PARAMETERS "ctcubemask::get_parameters()"
39 40
40 /* __ Debug definitions __________________________________________________ */ 41 /* __ Debug definitions __________________________________________________ */
41 42
@@ -336,6 +337,9 @@ void ctcubemask::save(void) @@ -336,6 +337,9 @@ void ctcubemask::save(void)
336 /***********************************************************************//** 337 /***********************************************************************//**
337 * @brief Get application parameters 338 * @brief Get application parameters
338 * 339 *
  340 + * @exception GException::invalid_value
  341 + * Parameter "inobs" is required for ctcubemask.
  342 + *
339 * Get all task parameters from parameter file or (if required) by querying 343 * Get all task parameters from parameter file or (if required) by querying
340 * the user. 344 * the user.
341 * 345 *
@@ -345,42 +349,23 @@ void ctcubemask::save(void) @@ -345,42 +349,23 @@ void ctcubemask::save(void)
345 ***************************************************************************/ 349 ***************************************************************************/
346 void ctcubemask::get_parameters(void) 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  
350 - if (m_obs.size() == 0) {  
351 -  
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); 352 + // If there are no observations in container then load them via user parameters
  353 + if (m_obs.size() == 0) {
363 354
364 - // Append CTA observation to container  
365 - m_obs.append(obs); 355 + // Throw exception if no infile is given, since this tool needs an observation
  356 + // including events
  357 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
366 358
367 - // Signal that no XML file should be used for storage  
368 - m_use_xml = false;  
369 -  
370 - }  
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 -  
378 - // Signal that XML file should be used for storage  
379 - m_use_xml = true; 359 + std::string msg = "Parameter \"inobs\" is required to be given in ctcubemask."
  360 + "Specify a vaild observation definition (XML or FITS) file to proceed"
  361 + "It is not possible to build a binned observation from scratch";
  362 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  363 + }
380 364
381 - } 365 + // Build observation container without response (not needed)
  366 + m_obs = get_observations(false);
382 367
383 - } // endif: there was no observation in the container 368 + } // endif: there was no observation in the container
384 369
385 // Get parameters 370 // Get parameters
386 m_regfile = (*this)["regfile"].filename(); 371 m_regfile = (*this)["regfile"].filename();
@@ -396,7 +381,7 @@ void ctcubemask::get_parameters(void) @@ -396,7 +381,7 @@ void ctcubemask::get_parameters(void)
396 // Optionally read ahead parameters so that they get correctly 381 // Optionally read ahead parameters so that they get correctly
397 // dumped into the log file 382 // dumped into the log file
398 if (read_ahead()) { 383 if (read_ahead()) {
399 - m_outfile = (*this)["outfile"].filename(); 384 + m_outcube = (*this)["outcube"].filename();
400 } 385 }
401 386
402 // Return 387 // Return
@@ -538,9 +523,9 @@ void ctcubemask::apply_mask(GCTAObservation* obs) @@ -538,9 +523,9 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
538 void ctcubemask::init_members(void) 523 void ctcubemask::init_members(void)
539 { 524 {
540 // Initialise parameters 525 // Initialise parameters
541 - m_infile.clear(); 526 + //m_infile.clear();
542 m_regfile.clear(); 527 m_regfile.clear();
543 - m_outfile.clear(); 528 + m_outcube.clear();
544 m_prefix.clear(); 529 m_prefix.clear();
545 m_usepnt = false; 530 m_usepnt = false;
546 m_ra = 0.0; 531 m_ra = 0.0;
@@ -552,7 +537,7 @@ void ctcubemask::init_members(void) @@ -552,7 +537,7 @@ void ctcubemask::init_members(void)
552 // Initialise protected members 537 // Initialise protected members
553 m_obs.clear(); 538 m_obs.clear();
554 m_infiles.clear(); 539 m_infiles.clear();
555 - m_use_xml = false; 540 +// m_use_xml = false;
556 541
557 // Return 542 // Return
558 return; 543 return;
@@ -567,9 +552,9 @@ void ctcubemask::init_members(void) @@ -567,9 +552,9 @@ void ctcubemask::init_members(void)
567 void ctcubemask::copy_members(const ctcubemask& app) 552 void ctcubemask::copy_members(const ctcubemask& app)
568 { 553 {
569 // Copy parameters 554 // Copy parameters
570 - m_infile = app.m_infile; 555 + //m_infile = app.m_infile;
571 m_regfile = app.m_regfile; 556 m_regfile = app.m_regfile;
572 - m_outfile = app.m_outfile; 557 + m_outcube = app.m_outcube;
573 m_prefix = app.m_prefix; 558 m_prefix = app.m_prefix;
574 m_usepnt = app.m_usepnt; 559 m_usepnt = app.m_usepnt;
575 m_ra = app.m_ra; 560 m_ra = app.m_ra;
@@ -581,7 +566,7 @@ void ctcubemask::copy_members(const ctcubemask&amp; app) @@ -581,7 +566,7 @@ void ctcubemask::copy_members(const ctcubemask&amp; app)
581 // Copy protected members 566 // Copy protected members
582 m_obs = app.m_obs; 567 m_obs = app.m_obs;
583 m_infiles = app.m_infiles; 568 m_infiles = app.m_infiles;
584 - m_use_xml = app.m_use_xml; 569 +// m_use_xml = app.m_use_xml;
585 570
586 // Return 571 // Return
587 return; 572 return;
@@ -622,7 +607,7 @@ std::string ctcubemask::check_infile(const std::string&amp; filename) const @@ -622,7 +607,7 @@ std::string ctcubemask::check_infile(const std::string&amp; filename) const
622 } 607 }
623 catch (GException::fits_hdu_not_found& e) { 608 catch (GException::fits_hdu_not_found& e) {
624 message = "No \"IMAGE\" extension found in input file \"" 609 message = "No \"IMAGE\" extension found in input file \""
625 - + m_outfile + "\"."; 610 + + m_outcube + "\".";
626 } 611 }
627 612
628 // Return 613 // Return
@@ -661,13 +646,13 @@ std::string ctcubemask::set_outfile_name(const std::string&amp; filename) const @@ -661,13 +646,13 @@ std::string ctcubemask::set_outfile_name(const std::string&amp; filename) const
661 void ctcubemask::save_fits(void) 646 void ctcubemask::save_fits(void)
662 { 647 {
663 // Get output filename 648 // Get output filename
664 - m_outfile = (*this)["outfile"].filename(); 649 + m_outcube = (*this)["outcube"].filename();
665 650
666 // Get CTA observation from observation container 651 // Get CTA observation from observation container
667 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]); 652 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
668 653
669 // Save event list 654 // Save event list
670 - save_counts_map(obs, m_outfile); 655 + save_counts_map(obs, m_outcube);
671 656
672 // Return 657 // Return
673 return; 658 return;
@@ -688,14 +673,14 @@ void ctcubemask::save_fits(void) @@ -688,14 +673,14 @@ void ctcubemask::save_fits(void)
688 void ctcubemask::save_xml(void) 673 void ctcubemask::save_xml(void)
689 { 674 {
690 // Get output filename and prefix 675 // Get output filename and prefix
691 - m_outfile = (*this)["outfile"].filename(); 676 + m_outcube = (*this)["m_outcube"].filename();
692 m_prefix = (*this)["prefix"].string(); 677 m_prefix = (*this)["prefix"].string();
693 678
694 // Issue warning if output filename has no .xml suffix 679 // Issue warning if output filename has no .xml suffix
695 - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); 680 + std::string suffix = gammalib::tolower(m_outcube.substr(m_outcube.length()-4,4));
696 if (suffix != ".xml") { 681 if (suffix != ".xml") {
697 log << "*** WARNING: Name of observation definition output file \""+ 682 log << "*** WARNING: Name of observation definition output file \""+
698 - m_outfile+"\"" << std::endl; 683 + m_outcube+"\"" << std::endl;
699 log << "*** WARNING: does not terminate with \".xml\"." << std::endl; 684 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
700 log << "*** WARNING: This is not an error, but might be misleading." 685 log << "*** WARNING: This is not an error, but might be misleading."
701 " It is recommended" << std::endl; 686 " It is recommended" << std::endl;
@@ -726,7 +711,7 @@ void ctcubemask::save_xml(void) @@ -726,7 +711,7 @@ void ctcubemask::save_xml(void)
726 } // endfor: looped over observations 711 } // endfor: looped over observations
727 712
728 // Save observations in XML file 713 // Save observations in XML file
729 - m_obs.save(m_outfile); 714 + m_obs.save(m_outcube);
730 715
731 // Return 716 // Return
732 return; 717 return;
src/ctcubemask/ctcubemask.hpp
@@ -36,7 +36,7 @@ @@ -36,7 +36,7 @@
36 36
37 /* __Definitions _________________________________________________________ */ 37 /* __Definitions _________________________________________________________ */
38 #define CTCUBEMASK_NAME "ctcubemask" 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,9 +78,9 @@ protected:
78 const std::string& outfile) const; 78 const std::string& outfile) const;
79 79
80 // User parameters 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 std::string m_regfile; //!< ds9 region file 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 std::string m_prefix; //!< Prefix for multiple counts maps 84 std::string m_prefix; //!< Prefix for multiple counts maps
85 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters 85 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
86 double m_ra; //!< RA of ROI centre 86 double m_ra; //!< RA of ROI centre
@@ -92,7 +92,7 @@ protected: @@ -92,7 +92,7 @@ protected:
92 // Protected members 92 // Protected members
93 GObservations m_obs; //!< Observations container 93 GObservations m_obs; //!< Observations container
94 std::vector<std::string> m_infiles; //!< Input event filenames 94 std::vector<std::string> m_infiles; //!< Input event filenames
95 - bool m_use_xml; //!< Use XML file instead of FITS file 95 + //bool m_use_xml; //!< Use XML file instead of FITS file
96 }; 96 };
97 97
98 98
src/ctcubemask/ctcubemask.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctcubemask.par # 7 # File: ctcubemask.par #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Chia-Chun Lu # 10 # Author: Chia-Chun Lu #
11 # c.lu@mpi-hd.mpg.de # 11 # c.lu@mpi-hd.mpg.de #
@@ -17,9 +17,9 @@ @@ -17,9 +17,9 @@
17 # 17 #
18 # File information 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 regfile, f, a, "NONE",,,"Exclusion region file in ds9 format" 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 prefix, s, h, "filtered_",,,"Prefix for counts cube in observation definition file" 23 prefix, s, h, "filtered_",,,"Prefix for counts cube in observation definition file"
24 # 24 #
25 # Selection parameters 25 # Selection parameters
src/ctexpcube/ctexpcube.cpp
@@ -284,10 +284,10 @@ void ctexpcube::save(void) @@ -284,10 +284,10 @@ void ctexpcube::save(void)
284 } 284 }
285 285
286 // Get output filename 286 // Get output filename
287 - m_outfile = (*this)["outfile"].filename(); 287 + m_outcube = (*this)["outcube"].filename();
288 288
289 // Save exposure cube 289 // Save exposure cube
290 - m_expcube.save(m_outfile, clobber()); 290 + m_expcube.save(m_outcube, clobber());
291 291
292 // Return 292 // Return
293 return; 293 return;
@@ -306,13 +306,13 @@ void ctexpcube::save(void) @@ -306,13 +306,13 @@ void ctexpcube::save(void)
306 void ctexpcube::init_members(void) 306 void ctexpcube::init_members(void)
307 { 307 {
308 // Initialise members 308 // Initialise members
309 - m_outfile.clear(); 309 + m_outcube.clear();
310 m_apply_edisp = false; 310 m_apply_edisp = false;
311 311
312 // Initialise protected members 312 // Initialise protected members
313 m_obs.clear(); 313 m_obs.clear();
314 m_expcube.clear(); 314 m_expcube.clear();
315 - m_ebounds.clear(); 315 + //m_ebounds.clear();
316 316
317 // Return 317 // Return
318 return; 318 return;
@@ -327,13 +327,12 @@ void ctexpcube::init_members(void) @@ -327,13 +327,12 @@ void ctexpcube::init_members(void)
327 void ctexpcube::copy_members(const ctexpcube& app) 327 void ctexpcube::copy_members(const ctexpcube& app)
328 { 328 {
329 // Copy attributes 329 // Copy attributes
330 - m_outfile = app.m_outfile; 330 + m_outcube = app.m_outcube;
331 m_apply_edisp = app.m_apply_edisp; 331 m_apply_edisp = app.m_apply_edisp;
332 332
333 // Copy protected members 333 // Copy protected members
334 m_obs = app.m_obs; 334 m_obs = app.m_obs;
335 m_expcube = app.m_expcube; 335 m_expcube = app.m_expcube;
336 - m_ebounds = app.m_ebounds;  
337 336
338 // Return 337 // Return
339 return; 338 return;
@@ -352,6 +351,8 @@ void ctexpcube::free_members(void) @@ -352,6 +351,8 @@ void ctexpcube::free_members(void)
352 351
353 /***********************************************************************//** 352 /***********************************************************************//**
354 * @brief Get application parameters 353 * @brief Get application parameters
  354 + * @exception GException::invalid_value
  355 + * Parameter "inobs" is required for ctexpcube.
355 * 356 *
356 * Get all task parameters from parameter file or (if required) by querying 357 * Get all task parameters from parameter file or (if required) by querying
357 * the user. The parameters are read in the correct order. 358 * the user. The parameters are read in the correct order.
@@ -360,159 +361,72 @@ void ctexpcube::free_members(void) @@ -360,159 +361,72 @@ void ctexpcube::free_members(void)
360 ***************************************************************************/ 361 ***************************************************************************/
361 void ctexpcube::get_parameters(void) 362 void ctexpcube::get_parameters(void)
362 { 363 {
363 - // If we do not have any observations in the container then get an  
364 - // input file name or observation descriptor file 364 + // If there are no observations in container then load them via user parameters
365 if (m_obs.size() == 0) { 365 if (m_obs.size() == 0) {
366 - get_obs();  
367 - }  
368 366
369 - // For all observations that have no response, set the response  
370 - // from the task parameters  
371 - set_response(m_obs); 367 + // Throw exception if no infile is given, since this tool needs an observation
  368 + // including events
  369 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
  370 +
  371 + std::string msg = "Parameter \"inobs\" is required to be given in ctexpcube."
  372 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  373 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  374 + }
  375 +
  376 + // Build observation container
  377 + m_obs = get_observations();
  378 +
  379 + } // endif: there was no observation in the container
372 380
373 // Read energy dispersion flag 381 // Read energy dispersion flag
374 m_apply_edisp = (*this)["edisp"].boolean(); 382 m_apply_edisp = (*this)["edisp"].boolean();
375 383
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(); 384 + // Get the incube filename
  385 + std::string incube = (*this)["incube"].filename();
392 386
393 - // Define exposure cube  
394 - m_expcube = GCTAExposure(wcs, coordsys, xref, yref,  
395 - -binsz, binsz, nxpix, nypix,  
396 - m_ebounds);  
397 - }  
398 -  
399 - // ... otherwise setup the exposure cube from the counts map  
400 - else { 387 + // Check for filename validity
  388 + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
401 389
402 - // Set exposure cube from counts map  
403 - set_from_cntmap(cntmap);  
404 -  
405 - } 390 + // Initialise event cube
  391 + GCTAEventCube cube;
406 392
407 - // Read output filename (if needed)  
408 - if (read_ahead()) {  
409 - m_outfile = (*this)["outfile"].filename();  
410 - } 393 + // Check for usepnt flag
  394 + bool usepnt = (*this)["usepnt"].boolean();
411 395
412 - // Return  
413 - return;  
414 -} 396 + // Check if pointing should be used
  397 + if (usepnt) {
415 398
  399 + // build cube from user parameters and pointing information
  400 + cube = get_cube(m_obs);
416 401
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 { 402 + } // endif: pointing used as map center
429 403
430 - // Allocate CTA observation  
431 - GCTAObservation obs;  
432 -  
433 - // Load input file in CTA observation  
434 - obs.load(filename); 404 + else {
435 405
436 - // Append CTA observation to container  
437 - m_obs.append(obs); 406 + // Build cube from user parameters
  407 + cube = get_cube();
438 408
439 -  
440 - }  
441 -  
442 - // ... otherwise try to open as XML file  
443 - catch (GException::fits_open_error &e) { 409 + } // endelse: m_usepnt was false
444 410
445 - // Load observations from XML file  
446 - m_obs.load(filename); 411 + // Define exposure cube
  412 + m_expcube = GCTAExposure(cube);
447 413
448 - } 414 + } // endif: filename was not valid
449 415
450 - // Return  
451 - return;  
452 -} 416 + // ... otherwise setup the exposure cube from the counts map
  417 + else {
453 418
  419 + // Set exposure cube from counts map
  420 + GCTAEventCube cube = set_from_cntmap(incube);
454 421
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 - } 422 + // Define exposure cube
  423 + m_expcube = GCTAExposure(cube);
508 424
509 - } // endif: observation contained an events cube 425 + } // endelse: cube loaded from file
510 426
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); 427 + // Read output filename (if needed)
  428 + if (read_ahead()) {
  429 + m_outcube = (*this)["outcube"].filename();
516 } 430 }
517 431
518 // Return 432 // Return
src/ctexpcube/ctexpcube.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTEXPCUBE_NAME "ctexpcube" 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,17 +67,14 @@ protected:
67 void copy_members(const ctexpcube& app); 67 void copy_members(const ctexpcube& app);
68 void free_members(void); 68 void free_members(void);
69 void get_parameters(void); 69 void get_parameters(void);
70 - void get_obs(void);  
71 - void set_from_cntmap(const std::string& filename);  
72 70
73 // User parameters 71 // User parameters
74 - std::string m_outfile; //!< Output exposure cube file 72 + std::string m_outcube; //!< Output exposure cube file
75 bool m_apply_edisp; //!< Apply energy dispersion? 73 bool m_apply_edisp; //!< Apply energy dispersion?
76 74
77 // Protected members 75 // Protected members
78 GObservations m_obs; //!< Observation container 76 GObservations m_obs; //!< Observation container
79 GCTAExposure m_expcube; //!< Exposure cube 77 GCTAExposure m_expcube; //!< Exposure cube
80 - GEbounds m_ebounds; //!< Energy boundaries  
81 }; 78 };
82 79
83 80
src/ctexpcube/ctexpcube.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctexpcube.par # 7 # File: ctexpcube.par #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Juergen Knoedlseder # 10 # Author: Juergen Knoedlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -15,11 +15,11 @@ @@ -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 caldb, s, a, "dummy",,, "Calibration database" 23 caldb, s, a, "dummy",,, "Calibration database"
24 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 edisp, b, h, no,,, "Apply energy dispersion?" 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,6 +34,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
34 # 34 #
35 # Spatial binning parameters 35 # Spatial binning parameters
36 #=========================== 36 #===========================
  37 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
37 nxpix, i, a, 200,,, "Size of the X axis in pixels" 38 nxpix, i, a, 200,,, "Size of the X axis in pixels"
38 nypix, i, a, 200,,, "Size of the Y axis in pixels" 39 nypix, i, a, 200,,, "Size of the Y axis in pixels"
39 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 40 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/ctlike/ctlike.cpp
@@ -326,11 +326,11 @@ void ctlike::save(void) @@ -326,11 +326,11 @@ void ctlike::save(void)
326 } 326 }
327 327
328 // Get output filename 328 // Get output filename
329 - m_outmdl = (*this)["outmdl"].filename(); 329 + m_outmodel = (*this)["outmodel"].filename();
330 330
331 // Write results out as XML model 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 // Return 336 // Return
@@ -341,6 +341,9 @@ void ctlike::save(void) @@ -341,6 +341,9 @@ void ctlike::save(void)
341 /***********************************************************************//** 341 /***********************************************************************//**
342 * @brief Get application parameters 342 * @brief Get application parameters
343 * 343 *
  344 + * @exception GException::invalid_value
  345 + * Parameter "inobs" is required for ctlike.
  346 + *
344 * Get all required task parameters from the parameter file or (if specified) 347 * Get all required task parameters from the parameter file or (if specified)
345 * by querying the user. Observation dependent parameters will only be read 348 * by querying the user. Observation dependent parameters will only be read
346 * if the observation container is actually empty. Observation dependent 349 * if the observation container is actually empty. Observation dependent
@@ -357,56 +360,38 @@ void ctlike::save(void) @@ -357,56 +360,38 @@ void ctlike::save(void)
357 ***************************************************************************/ 360 ***************************************************************************/
358 void ctlike::get_parameters(void) 361 void ctlike::get_parameters(void)
359 { 362 {
360 - // If there are no observations in container then add a single CTA  
361 - // observation using the parameters from the parameter file 363 + // If there are no observations in container then load them via user parameters
362 if (m_obs.size() == 0) { 364 if (m_obs.size() == 0) {
363 365
364 - // Allocate CTA observation  
365 - GCTAObservation obs;  
366 -  
367 - // Get event file name  
368 - std::string filename = (*this)["infile"].filename();  
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 - 366 + // Throw exception if no infile is given, since this tool needs an observation
  367 + // including events
  368 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
  369 + std::string msg = "Parameter \"inobs\" is required to be given in ctlike."
  370 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  371 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
388 } 372 }
389 373
390 - // ... otherwise try to open as XML file  
391 - catch (GException::fits_open_error &e) { 374 + // Build observation container
  375 + m_obs = get_observations();
392 376
393 - // Load observations from XML file  
394 - m_obs.load(filename); 377 + } // endif: there was no observation in the container
395 378
396 - // For all observations that have no response, set the response  
397 - // from the task parameters  
398 - set_response(m_obs); 379 + // If only single observation is used, read statistics parameter
  380 + if (!m_use_xml) {
399 381
400 - } // endcatch: file was an XML file 382 + // Get other task parameters
  383 + std::string statistics = gammalib::toupper((*this)["stat"].string());
401 384
402 - } // endif: there was no observation in the container 385 + // Set statistics
  386 + (*m_obs[0]).statistics(statistics);
  387 + }
403 388
404 // If there is are no models associated with the observations then 389 // If there is are no models associated with the observations then
405 // load now the model definition 390 // load now the model definition
406 if (m_obs.models().size() == 0) { 391 if (m_obs.models().size() == 0) {
407 392
408 // Get models XML filename 393 // Get models XML filename
409 - std::string filename = (*this)["srcmdl"].filename(); 394 + std::string filename = (*this)["inmodel"].filename();
410 395
411 // Setup models for optimizing. 396 // Setup models for optimizing.
412 m_obs.models(GModels(filename)); 397 m_obs.models(GModels(filename));
@@ -420,7 +405,7 @@ void ctlike::get_parameters(void) @@ -420,7 +405,7 @@ void ctlike::get_parameters(void)
420 // Optionally read ahead parameters so that they get correctly 405 // Optionally read ahead parameters so that they get correctly
421 // dumped into the log file 406 // dumped into the log file
422 if (read_ahead()) { 407 if (read_ahead()) {
423 - m_outmdl = (*this)["outmdl"].filename(); 408 + m_outmodel = (*this)["outmodel"].filename();
424 } 409 }
425 410
426 // Return 411 // Return
@@ -547,7 +532,7 @@ double ctlike::reoptimize_lm(void) @@ -547,7 +532,7 @@ double ctlike::reoptimize_lm(void)
547 void ctlike::init_members(void) 532 void ctlike::init_members(void)
548 { 533 {
549 // Initialise members 534 // Initialise members
550 - m_outmdl.clear(); 535 + m_outmodel.clear();
551 m_obs.clear(); 536 m_obs.clear();
552 m_refit = false; 537 m_refit = false;
553 m_max_iter = 100; // Set maximum number of iterations 538 m_max_iter = 100; // Set maximum number of iterations
@@ -573,7 +558,7 @@ void ctlike::copy_members(const ctlike&amp; app) @@ -573,7 +558,7 @@ void ctlike::copy_members(const ctlike&amp; app)
573 { 558 {
574 // Copy attributes 559 // Copy attributes
575 m_refit = app.m_refit; 560 m_refit = app.m_refit;
576 - m_outmdl = app.m_outmdl; 561 + m_outmodel = app.m_outmodel;
577 m_obs = app.m_obs; 562 m_obs = app.m_obs;
578 m_max_iter = app.m_max_iter; 563 m_max_iter = app.m_max_iter;
579 m_max_stall = app.m_max_stall; 564 m_max_stall = app.m_max_stall;
src/ctlike/ctlike.hpp
@@ -35,7 +35,7 @@ @@ -35,7 +35,7 @@
35 35
36 /* __Definitions _________________________________________________________ */ 36 /* __Definitions _________________________________________________________ */
37 #define CTLIKE_NAME "ctlike" 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,7 +74,7 @@ protected:
74 74
75 // User parameters 75 // User parameters
76 bool m_refit; //!< Refitting 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 bool m_apply_edisp; //!< Apply energy dispersion? 78 bool m_apply_edisp; //!< Apply energy dispersion?
79 79
80 // Members 80 // Members
src/ctlike/ctlike.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctlike.par # 7 # File: ctlike.par #
8 -# Version: 00-05-00 # 8 +# Version: 00-06-00 #
9 # # 9 # #
10 # Author: Jurgen Knodlseder # 10 # Author: Jurgen Knodlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -17,14 +17,14 @@ @@ -17,14 +17,14 @@
17 # 17 #
18 # General parameters 18 # General parameters
19 #=================== 19 #===================
20 -infile, f, a, "events.fits",,, "Event list, counts cube or observation definition file" 20 +inobs, f, a, "events.fits",,, "Event list, counts cube or observation definition file"
21 stat, s, h, "POISSON",,, "Optimization statistics" 21 stat, s, h, "POISSON",,, "Optimization statistics"
22 edisp, b, h, no,,, "Apply energy dispersion?" 22 edisp, b, h, no,,, "Apply energy dispersion?"
23 refit, b, h, no,,, "Do refitting?" 23 refit, b, h, no,,, "Do refitting?"
24 caldb, s, a, "dummy",,, "Calibration database" 24 caldb, s, a, "dummy",,, "Calibration database"
25 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 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" 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 # Standard parameters 29 # Standard parameters
30 #==================== 30 #====================
src/ctmodel/ctmodel.cpp
@@ -207,9 +207,6 @@ void ctmodel::run(void) @@ -207,9 +207,6 @@ void ctmodel::run(void)
207 // Get task parameters 207 // Get task parameters
208 get_parameters(); 208 get_parameters();
209 209
210 - // Setup observation container  
211 - setup_obs();  
212 -  
213 // Write parameters into logger 210 // Write parameters into logger
214 if (logTerse()) { 211 if (logTerse()) {
215 log_parameters(); 212 log_parameters();
@@ -246,11 +243,6 @@ void ctmodel::run(void) @@ -246,11 +243,6 @@ void ctmodel::run(void)
246 log.header1("Generate model cube"); 243 log.header1("Generate model cube");
247 } 244 }
248 245
249 - // Initialise counts cube if cube has not yet been set  
250 - if (!has_cube()) {  
251 - init_cube();  
252 - }  
253 -  
254 // Loop over all observations in the container 246 // Loop over all observations in the container
255 for (int i = 0; i < m_obs.size(); ++i) { 247 for (int i = 0; i < m_obs.size(); ++i) {
256 248
@@ -318,10 +310,10 @@ void ctmodel::save(void) @@ -318,10 +310,10 @@ void ctmodel::save(void)
318 } 310 }
319 311
320 // Make sure we have the FITS filename 312 // Make sure we have the FITS filename
321 - m_outfile = (*this)["outfile"].filename(); 313 + m_outcube = (*this)["outcube"].filename();
322 314
323 // Save model cube into FITS file 315 // Save model cube into FITS file
324 - m_cube.save(m_outfile, clobber()); 316 + m_cube.save(m_outcube, clobber());
325 317
326 // Return 318 // Return
327 return; 319 return;
@@ -365,35 +357,12 @@ void ctmodel::cube(const GCTAEventCube&amp; cube) @@ -365,35 +357,12 @@ void ctmodel::cube(const GCTAEventCube&amp; cube)
365 void ctmodel::init_members(void) 357 void ctmodel::init_members(void)
366 { 358 {
367 // Initialise members 359 // 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; 360 + m_outcube.clear();
391 m_apply_edisp = false; 361 m_apply_edisp = false;
392 362
393 // Initialise protected members 363 // Initialise protected members
394 m_obs.clear(); 364 m_obs.clear();
395 m_cube.clear(); 365 m_cube.clear();
396 - m_ebounds.clear();  
397 m_gti.clear(); 366 m_gti.clear();
398 m_has_cube = false; 367 m_has_cube = false;
399 368
@@ -410,35 +379,12 @@ void ctmodel::init_members(void) @@ -410,35 +379,12 @@ void ctmodel::init_members(void)
410 void ctmodel::copy_members(const ctmodel& app) 379 void ctmodel::copy_members(const ctmodel& app)
411 { 380 {
412 // Copy attributes 381 // 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; 382 + m_outcube = app.m_outcube;
436 m_apply_edisp = app.m_apply_edisp; 383 m_apply_edisp = app.m_apply_edisp;
437 384
438 // Copy protected members 385 // Copy protected members
439 m_obs = app.m_obs; 386 m_obs = app.m_obs;
440 m_cube = app.m_cube; 387 m_cube = app.m_cube;
441 - m_ebounds = app.m_ebounds;  
442 m_gti = app.m_gti; 388 m_gti = app.m_gti;
443 m_has_cube = app.m_has_cube; 389 m_has_cube = app.m_has_cube;
444 390
@@ -465,255 +411,83 @@ void ctmodel::free_members(void) @@ -465,255 +411,83 @@ void ctmodel::free_members(void)
465 ***************************************************************************/ 411 ***************************************************************************/
466 void ctmodel::get_parameters(void) 412 void ctmodel::get_parameters(void)
467 { 413 {
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();  
474 - } 414 + // If there are no observations in container then load them via user parameters
  415 + if (m_obs.size() == 0) {
475 416
476 - // Read model definiton file if required  
477 - if (m_obs.models().size() == 0) {  
478 - m_srcmdl = (*this)["srcmdl"].filename();  
479 - } 417 + // Build observation container
  418 + m_obs = get_observations();
480 419
481 - // Get energy dispersion flag parameters  
482 - m_apply_edisp = (*this)["edisp"].boolean(); 420 + } // endif: there was no observation in the container
483 421
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 - } 422 + // Read model definition file if required
  423 + if (m_obs.models().size() == 0) {
  424 + std::string inmodel = (*this)["inmodel"].filename();
503 425
504 - } 426 + // Load models from file
  427 + m_obs.models(inmodel);
505 428
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 - } 429 + } // endif: there were no models
517 430
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 - }  
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(); 431 + // Read optionally output cube filenames
  432 + if (read_ahead()) {
  433 + m_outcube = (*this)["outcube"].filename();
538 } 434 }
539 435
540 - // Return  
541 - return;  
542 -}  
543 -  
544 -  
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  
616 - else {  
617 -  
618 - // First try to open the file as an event list or counts map  
619 - try { 436 + // Get energy dispersion flag parameters
  437 + m_apply_edisp = (*this)["edisp"].boolean();
620 438
621 - // Allocate CTA observation  
622 - GCTAObservation obs; 439 + if (!m_has_cube) {
623 440
624 - // Load CTA observation from FITS file  
625 - obs.load(m_obsfile); 441 + // Read cube definition file
  442 + std::string incube = (*this)["incube"].filename();
626 443
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) { 444 + // If no cube file has been specified then read all parameters that
  445 + // are necessary to create the cube from scratch (see method ctool::get_cube())
  446 + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
634 447
635 - // Load observations from XML file. This will throw  
636 - // an exception if it fails.  
637 - m_obs.load(m_obsfile); 448 + // Get parameters
  449 + bool usepnt = (*this)["usepnt"].boolean();
638 450
639 - } 451 + // Initialise event cube
  452 + GCTAEventCube cube;
640 453
641 - } // endelse: loaded information from input file 454 + // Check if pointing should be used
  455 + if (usepnt) {
642 456
643 - } // endif: there was no observation in the container 457 + // buld cube from user parameters and pointing information
  458 + m_cube = get_cube(m_obs);
644 459
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 - } 460 + } // endif: pointing used as map center
650 461
651 - // For all observations that have no response, set the response  
652 - // from the task parameters  
653 - set_response(m_obs); 462 + else {
654 463
655 - // Return  
656 - return;  
657 -} 464 + // Build cube from user parameters
  465 + m_cube = get_cube();
658 466
  467 + } // endelse: m_usepnt was false
659 468
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 - } 469 + } // endif: incube file was not vaild
695 470
696 - // ... otherwise load cube from file  
697 - else { 471 + // ... else: build map information from file
  472 + else {
698 473
699 - // Load cube  
700 - m_cube.load(m_infile); 474 + // Load cube from given file
  475 + m_cube.load(incube);
701 476
702 - // Set all cube bins to zero  
703 - for (int i = 0; i < m_cube.size(); ++i) {  
704 - m_cube[i]->counts(0.0); 477 + // Set all cube bins to zero
  478 + for (int i = 0; i < m_cube.size(); ++i) {
  479 + m_cube[i]->counts(0.0);
  480 + }
705 } 481 }
706 482
  483 + // Signal that cube has been set
  484 + m_has_cube = true;
707 } 485 }
708 486
709 - // Signal that cube has been set  
710 - m_has_cube = true;  
711 -  
712 // Return 487 // Return
713 return; 488 return;
714 } 489 }
715 490
716 -  
717 /***********************************************************************//** 491 /***********************************************************************//**
718 * @brief Fill model into model cube 492 * @brief Fill model into model cube
719 * 493 *
src/ctmodel/ctmodel.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTMODEL_NAME "ctmodel" 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,41 +78,16 @@ protected:
78 void copy_members(const ctmodel& app); 78 void copy_members(const ctmodel& app);
79 void free_members(void); 79 void free_members(void);
80 void get_parameters(void); 80 void get_parameters(void);
81 - void setup_obs(void);  
82 - void init_cube(void);  
83 void fill_cube(const GCTAObservation* obs); 81 void fill_cube(const GCTAObservation* obs);
84 bool has_cube(void) const; 82 bool has_cube(void) const;
85 83
86 // User parameters 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 bool m_apply_edisp; //!< Apply energy dispersion? 86 bool m_apply_edisp; //!< Apply energy dispersion?
111 87
112 // Protected members 88 // Protected members
113 GObservations m_obs; //!< Observation container 89 GObservations m_obs; //!< Observation container
114 GCTAEventCube m_cube; //!< Model cube 90 GCTAEventCube m_cube; //!< Model cube
115 - GEbounds m_ebounds; //!< Energy boundaries  
116 GGti m_gti; //!< Model cube GTIs 91 GGti m_gti; //!< Model cube GTIs
117 bool m_has_cube; //!< Signal if cube has been set or loaded 92 bool m_has_cube; //!< Signal if cube has been set or loaded
118 }; 93 };
src/ctmodel/ctmodel.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctmodel.par # 7 # File: ctmodel.par #
8 -# Version: 00-03-00 # 8 +# Version: 00-04-00 #
9 # # 9 # #
10 # Author: Juergen Knoedlseder # 10 # Author: Juergen Knoedlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -17,18 +17,19 @@ @@ -17,18 +17,19 @@
17 # 17 #
18 # General parameters 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 edisp, b, h, no,,, "Apply energy dispersion?" 23 edisp, b, h, no,,, "Apply energy dispersion?"
24 # 24 #
25 # Observation parameters 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 caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" 28 caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
29 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 29 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
30 ra, r, a, 83.63,0,360, "RA of pointing (degrees)" 30 ra, r, a, 83.63,0,360, "RA of pointing (degrees)"
31 dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)" 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 tmin, r, a, 0.0,0,, "Start time (MET in s)" 33 tmin, r, a, 0.0,0,, "Start time (MET in s)"
33 tmax, r, a, 1800.0,0,, "End time (MET in s)" 34 tmax, r, a, 1800.0,0,, "End time (MET in s)"
34 deadc, r, h, 0.95,0,1, "Average deadtime correction factor" 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,6 +44,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
43 # 44 #
44 # Spatial binning parameters 45 # Spatial binning parameters
45 #=========================== 46 #===========================
  47 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
46 nxpix, i, a, 200,,, "Size of the X axis in pixels" 48 nxpix, i, a, 200,,, "Size of the X axis in pixels"
47 nypix, i, a, 200,,, "Size of the Y axis in pixels" 49 nypix, i, a, 200,,, "Size of the Y axis in pixels"
48 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 50 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/ctobssim/ctobssim.cpp
@@ -377,7 +377,7 @@ void ctobssim::run(void) @@ -377,7 +377,7 @@ void ctobssim::run(void)
377 outfile = m_prefix + gammalib::str(i) + ".fits"; 377 outfile = m_prefix + gammalib::str(i) + ".fits";
378 } 378 }
379 else { 379 else {
380 - outfile = (*this)["outfile"].filename(); 380 + outfile = (*this)["outevents"].filename();
381 } 381 }
382 382
383 // Store output file name in original observation 383 // Store output file name in original observation
@@ -427,7 +427,6 @@ void ctobssim::run(void) @@ -427,7 +427,6 @@ void ctobssim::run(void)
427 return; 427 return;
428 } 428 }
429 429
430 -  
431 /***********************************************************************//** 430 /***********************************************************************//**
432 * @brief Save the selected event list(s) 431 * @brief Save the selected event list(s)
433 * 432 *
@@ -490,35 +489,20 @@ void ctobssim::save(void) @@ -490,35 +489,20 @@ void ctobssim::save(void)
490 void ctobssim::init_members(void) 489 void ctobssim::init_members(void)
491 { 490 {
492 // Initialise user parameters 491 // Initialise user parameters
493 - m_infile.clear();  
494 - m_outfile.clear(); 492 + m_outevents.clear();
495 m_prefix.clear(); 493 m_prefix.clear();
496 m_seed = 1; 494 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;  
504 m_apply_edisp = false; 495 m_apply_edisp = false;
505 - m_deadc = 1.0;  
506 496
507 // Initialise protected members 497 // Initialise protected members
508 m_rans.clear(); 498 m_rans.clear();
509 m_obs.clear(); 499 m_obs.clear();
510 - m_use_xml = false;  
511 m_save_and_dispose = false; 500 m_save_and_dispose = false;
512 501
513 // Set fixed parameters 502 // Set fixed parameters
514 m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2 503 m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2
515 m_max_photons = 1000000; //!< Maximum number of photons / time slice 504 m_max_photons = 1000000; //!< Maximum number of photons / time slice
516 505
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 // Initialise first event identifier 506 // Initialise first event identifier
523 m_event_id = 1; 507 m_event_id = 1;
524 508
@@ -535,28 +519,17 @@ void ctobssim::init_members(void) @@ -535,28 +519,17 @@ void ctobssim::init_members(void)
535 void ctobssim::copy_members(const ctobssim& app) 519 void ctobssim::copy_members(const ctobssim& app)
536 { 520 {
537 // Copy user parameters 521 // Copy user parameters
538 - m_infile = app.m_infile;  
539 - m_outfile = app.m_outfile; 522 + m_outevents = app.m_outevents;
540 m_prefix = app.m_prefix; 523 m_prefix = app.m_prefix;
541 m_seed = app.m_seed; 524 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 m_apply_edisp = app.m_apply_edisp; 525 m_apply_edisp = app.m_apply_edisp;
550 - m_deadc = app.m_deadc;  
551 526
552 // Copy protected members 527 // Copy protected members
553 m_area = app.m_area; 528 m_area = app.m_area;
554 m_max_photons = app.m_max_photons; 529 m_max_photons = app.m_max_photons;
555 m_rans = app.m_rans; 530 m_rans = app.m_rans;
556 m_obs = app.m_obs; 531 m_obs = app.m_obs;
557 - m_use_xml = app.m_use_xml;  
558 m_save_and_dispose = app.m_save_and_dispose; 532 m_save_and_dispose = app.m_save_and_dispose;
559 - m_cta_ref = app.m_cta_ref;  
560 m_event_id = app.m_event_id; 533 m_event_id = app.m_event_id;
561 534
562 // Return 535 // Return
@@ -585,48 +558,22 @@ void ctobssim::free_members(void) @@ -585,48 +558,22 @@ void ctobssim::free_members(void)
585 ***************************************************************************/ 558 ***************************************************************************/
586 void ctobssim::get_parameters(void) 559 void ctobssim::get_parameters(void)
587 { 560 {
588 - // If there are no observations in container then add a single CTA  
589 - // observation using the parameters from the parameter file  
590 - if (m_obs.size() == 0) {  
591 -  
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); 561 + // If there are no observations in container then load them via user parameters
  562 + if (m_obs.size() == 0) {
614 563
615 - // Append CTA observation to container  
616 - m_obs.append(obs); 564 + // Build observation container
  565 + m_obs = get_observations();
617 566
618 - // Load models into container  
619 - m_obs.models(m_infile); 567 + } // endif: there was no observation in the container
620 568
621 - // Signal that no XML file should be used for storage  
622 - m_use_xml = false; 569 + // Read model definition file if required
  570 + if (m_obs.models().size() == 0) {
  571 + std::string inmodel = (*this)["inmodel"].filename();
623 572
624 - } // endif: there was no observation in the container 573 + // Load models from file
  574 + m_obs.models(inmodel);
625 575
626 - // ... otherwise we signal that an XML file should be written  
627 - else {  
628 - m_use_xml = true;  
629 - } 576 + } // endif: there were no models
630 577
631 // Get other parameters 578 // Get other parameters
632 m_seed = (*this)["seed"].integer(); 579 m_seed = (*this)["seed"].integer();
@@ -635,20 +582,10 @@ void ctobssim::get_parameters(void) @@ -635,20 +582,10 @@ void ctobssim::get_parameters(void)
635 // Optionally read ahead parameters so that they get correctly 582 // Optionally read ahead parameters so that they get correctly
636 // dumped into the log file 583 // dumped into the log file
637 if (read_ahead()) { 584 if (read_ahead()) {
638 - m_outfile = (*this)["outfile"].filename(); 585 + m_outevents = (*this)["outevents"].filename();
639 m_prefix = (*this)["prefix"].string(); 586 m_prefix = (*this)["prefix"].string();
640 } 587 }
641 588
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 - }  
650 - }  
651 -  
652 // Initialise random number generators. We initialise here one random 589 // Initialise random number generators. We initialise here one random
653 // number generator per observation so that each observation will 590 // number generator per observation so that each observation will
654 // get it's own random number generator. This will lead to identical 591 // get it's own random number generator. This will lead to identical
@@ -700,81 +637,6 @@ void ctobssim::get_parameters(void) @@ -700,81 +637,6 @@ void ctobssim::get_parameters(void)
700 637
701 638
702 /***********************************************************************//** 639 /***********************************************************************//**
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 * @brief Simulate source events from photon list 640 * @brief Simulate source events from photon list
779 * 641 *
780 * @param[in] obs Pointer on CTA observation. 642 * @param[in] obs Pointer on CTA observation.
@@ -1163,13 +1025,13 @@ void ctobssim::save_fits(void) @@ -1163,13 +1025,13 @@ void ctobssim::save_fits(void)
1163 if (!m_save_and_dispose) { 1025 if (!m_save_and_dispose) {
1164 1026
1165 // Get output filename 1027 // Get output filename
1166 - m_outfile = (*this)["outfile"].filename(); 1028 + m_outevents = (*this)["outevents"].filename();
1167 1029
1168 // Get CTA observation from observation container 1030 // Get CTA observation from observation container
1169 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]); 1031 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1170 1032
1171 // Save observation into FITS file 1033 // Save observation into FITS file
1172 - obs->save(m_outfile, clobber()); 1034 + obs->save(m_outevents, clobber());
1173 1035
1174 } // endif: event list has not yet been saved and disposed 1036 } // endif: event list has not yet been saved and disposed
1175 1037
@@ -1191,14 +1053,14 @@ void ctobssim::save_fits(void) @@ -1191,14 +1053,14 @@ void ctobssim::save_fits(void)
1191 void ctobssim::save_xml(void) 1053 void ctobssim::save_xml(void)
1192 { 1054 {
1193 // Get output filename and prefix 1055 // Get output filename and prefix
1194 - m_outfile = (*this)["outfile"].filename(); 1056 + m_outevents = (*this)["outevents"].filename();
1195 m_prefix = (*this)["prefix"].string(); 1057 m_prefix = (*this)["prefix"].string();
1196 1058
1197 // Issue warning if output filename has no .xml suffix 1059 // Issue warning if output filename has no .xml suffix
1198 - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); 1060 + std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4));
1199 if (suffix != ".xml") { 1061 if (suffix != ".xml") {
1200 log << "*** WARNING: Name of observation definition output file \""+ 1062 log << "*** WARNING: Name of observation definition output file \""+
1201 - m_outfile+"\"" << std::endl; 1063 + m_outevents+"\"" << std::endl;
1202 log << "*** WARNING: does not terminate with \".xml\"." << std::endl; 1064 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1203 log << "*** WARNING: This is not an error, but might be misleading." 1065 log << "*** WARNING: This is not an error, but might be misleading."
1204 " It is recommended" << std::endl; 1066 " It is recommended" << std::endl;
@@ -1239,7 +1101,7 @@ void ctobssim::save_xml(void) @@ -1239,7 +1101,7 @@ void ctobssim::save_xml(void)
1239 } // endif: event list has not yet been saved and disposed 1101 } // endif: event list has not yet been saved and disposed
1240 1102
1241 // Save observations in XML file 1103 // Save observations in XML file
1242 - m_obs.save(m_outfile); 1104 + m_obs.save(m_outevents);
1243 1105
1244 // Return 1106 // Return
1245 return; 1107 return;
src/ctobssim/ctobssim.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTOBSSIM_NAME "ctobssim" 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,7 +73,6 @@ protected:
73 void copy_members(const ctobssim& app); 73 void copy_members(const ctobssim& app);
74 void free_members(void); 74 void free_members(void);
75 void get_parameters(void); 75 void get_parameters(void);
76 - void set_list(GCTAObservation* obs);  
77 void simulate_source(GCTAObservation* obs, 76 void simulate_source(GCTAObservation* obs,
78 const GModels& models, 77 const GModels& models,
79 GRan& ran, 78 GRan& ran,
@@ -86,19 +85,10 @@ protected: @@ -86,19 +85,10 @@ protected:
86 void save_xml(void); 85 void save_xml(void);
87 86
88 // User parameters 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 std::string m_prefix; //!< Prefix for multiple event lists 89 std::string m_prefix; //!< Prefix for multiple event lists
92 int m_seed; //!< Random number generator seed 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 bool m_apply_edisp; //!< Apply energy dispersion? 91 bool m_apply_edisp; //!< Apply energy dispersion?
101 - double m_deadc; //!< Average deadtime correction  
102 92
103 // Protected members 93 // Protected members
104 mutable GObservations m_obs; //!< Observation container 94 mutable GObservations m_obs; //!< Observation container
@@ -106,8 +96,6 @@ protected: @@ -106,8 +96,6 @@ protected:
106 double m_area; //!< Surface area for simulation (cm2) 96 double m_area; //!< Surface area for simulation (cm2)
107 int m_max_photons; //!< Maximum number of photons/slice 97 int m_max_photons; //!< Maximum number of photons/slice
108 std::vector<GRan> m_rans; //!< Random number generators 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 int m_event_id; //!< Event identifier 99 int m_event_id; //!< Event identifier
112 }; 100 };
113 101
@@ -120,7 +108,6 @@ protected: @@ -120,7 +108,6 @@ protected:
120 inline 108 inline
121 const GObservations& ctobssim::obs(void) const 109 const GObservations& ctobssim::obs(void) const
122 { 110 {
123 - // Return observation container  
124 return m_obs; 111 return m_obs;
125 } 112 }
126 113
src/ctobssim/ctobssim.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctobssim.par # 7 # File: ctobssim.par #
8 -# Version: 00-08-00 # 8 +# Version: 00-10-00 #
9 # # 9 # #
10 # Author: Jurgen Knodlseder # 10 # Author: Jurgen Knodlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -17,8 +17,9 @@ @@ -17,8 +17,9 @@
17 # 17 #
18 # File information 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" 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"
22 prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" 23 prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file"
23 # 24 #
24 # Observation parameters 25 # Observation parameters
src/ctpsfcube/ctpsfcube.cpp
@@ -275,10 +275,10 @@ void ctpsfcube::save(void) @@ -275,10 +275,10 @@ void ctpsfcube::save(void)
275 } 275 }
276 276
277 // Get output filename 277 // Get output filename
278 - m_outfile = (*this)["outfile"].filename(); 278 + m_outcube = (*this)["outcube"].filename();
279 279
280 // Save PSF cube 280 // Save PSF cube
281 - m_psfcube.save(m_outfile, clobber()); 281 + m_psfcube.save(m_outcube, clobber());
282 282
283 // Return 283 // Return
284 return; 284 return;
@@ -297,13 +297,12 @@ void ctpsfcube::save(void) @@ -297,13 +297,12 @@ void ctpsfcube::save(void)
297 void ctpsfcube::init_members(void) 297 void ctpsfcube::init_members(void)
298 { 298 {
299 // Initialise members 299 // Initialise members
300 - m_outfile.clear(); 300 + m_outcube.clear();
301 m_apply_edisp = false; 301 m_apply_edisp = false;
302 302
303 // Initialise protected members 303 // Initialise protected members
304 m_obs.clear(); 304 m_obs.clear();
305 m_psfcube.clear(); 305 m_psfcube.clear();
306 - m_ebounds.clear();  
307 306
308 // Return 307 // Return
309 return; 308 return;
@@ -318,13 +317,12 @@ void ctpsfcube::init_members(void) @@ -318,13 +317,12 @@ void ctpsfcube::init_members(void)
318 void ctpsfcube::copy_members(const ctpsfcube& app) 317 void ctpsfcube::copy_members(const ctpsfcube& app)
319 { 318 {
320 // Copy attributes 319 // Copy attributes
321 - m_outfile = app.m_outfile; 320 + m_outcube = app.m_outcube;
322 m_apply_edisp = app.m_apply_edisp; 321 m_apply_edisp = app.m_apply_edisp;
323 322
324 // Copy protected members 323 // Copy protected members
325 m_obs = app.m_obs; 324 m_obs = app.m_obs;
326 m_psfcube = app.m_psfcube; 325 m_psfcube = app.m_psfcube;
327 - m_ebounds = app.m_ebounds;  
328 326
329 // Return 327 // Return
330 return; 328 return;
@@ -344,6 +342,9 @@ void ctpsfcube::free_members(void) @@ -344,6 +342,9 @@ void ctpsfcube::free_members(void)
344 /***********************************************************************//** 342 /***********************************************************************//**
345 * @brief Get application parameters 343 * @brief Get application parameters
346 * 344 *
  345 + * @exception GException::invalid_value
  346 + * Parameter "inobs" is required for ctbin.
  347 + *
347 * Get all task parameters from parameter file or (if required) by querying 348 * Get all task parameters from parameter file or (if required) by querying
348 * the user. The parameters are read in the correct order. 349 * the user. The parameters are read in the correct order.
349 * 350 *
@@ -351,162 +352,75 @@ void ctpsfcube::free_members(void) @@ -351,162 +352,75 @@ void ctpsfcube::free_members(void)
351 ***************************************************************************/ 352 ***************************************************************************/
352 void ctpsfcube::get_parameters(void) 353 void ctpsfcube::get_parameters(void)
353 { 354 {
354 - // If we do not have any observations in the container then get an  
355 - // input file name or observation descriptor file 355 + // If there are no observations in container then load them via user parameters
356 if (m_obs.size() == 0) { 356 if (m_obs.size() == 0) {
357 - get_obs();  
358 - }  
359 357
360 - // Make sure that response is set  
361 - set_response(m_obs); 358 + // Throw exception if no infile is given, since this tool needs an observation
  359 + // including events
  360 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
362 361
363 - // Read energy dispersion flag  
364 - m_apply_edisp = (*this)["edisp"].boolean(); 362 + std::string msg = "Parameter \"inobs\" is required to be given in ctpsfcube."
  363 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  364 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  365 + }
365 366
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 - } 367 + // Build observation container
  368 + m_obs = get_observations();
399 369
400 - // Read output filename (if needed)  
401 - if (read_ahead()) {  
402 - m_outfile = (*this)["outfile"].filename();  
403 - } 370 + } // endif: there was no observation in the container
404 371
405 - // Return  
406 - return;  
407 -} 372 + // Get the incube filename
  373 + std::string incube = (*this)["incube"].filename();
408 374
  375 + // Get additional binning parameters
  376 + double dmax = (*this)["amax"].real();
  377 + int ndbins = (*this)["anumbins"].integer();
409 378
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(); 379 + // Check for filename validity
  380 + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
419 381
420 - // Try first to open as FITS file  
421 - try { 382 + // Get parameters
  383 + bool usepnt = (*this)["usepnt"].boolean();
422 384
423 - // Allocate CTA observation  
424 - GCTAObservation obs;  
425 -  
426 - // Load input file in CTA observation  
427 - obs.load(filename); 385 + // Initialise event cube
  386 + GCTAEventCube cube;
428 387
429 - // Append CTA observation to container  
430 - m_obs.append(obs); 388 + // Check if pointing should be used
  389 + if (usepnt) {
431 390
432 - }  
433 -  
434 - // ... otherwise try to open as XML file  
435 - catch (GException::fits_open_error &e) { 391 + // buld cube from user parameters and pointing information
  392 + cube = get_cube(m_obs);
436 393
437 - // Load observations from XML file  
438 - m_obs.load(filename); 394 + } // endif: pointing used as map center
  395 +
  396 + else {
439 397
  398 + // Build cube from user parameters
  399 + cube = get_cube();
  400 +
  401 + } // endelse: m_usepnt was false
  402 +
  403 + // Define psf cube
  404 + m_psfcube = GCTAMeanPsf(cube, dmax, ndbins);
440 } 405 }
441 406
442 - // Return  
443 - return;  
444 -} 407 + // ... otherwise setup the exposure cube from the counts map
  408 + else {
445 409
  410 + // Set exposure cube from counts map
  411 + GCTAEventCube cube = set_from_cntmap(incube);
446 412
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 - } 413 + // Define psf cube
  414 + m_psfcube = GCTAMeanPsf(cube, dmax, ndbins);
502 415
503 - } // endif: observation contained an events cube 416 + } // endelse: cube loaded from file
504 417
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); 418 + // Read energy dispersion flag
  419 + m_apply_edisp = (*this)["edisp"].boolean();
  420 +
  421 + // Read output filename (if needed)
  422 + if (read_ahead()) {
  423 + m_outcube = (*this)["outcube"].filename();
510 } 424 }
511 425
512 // Return 426 // Return
src/ctpsfcube/ctpsfcube.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTPSFCUBE_NAME "ctpsfcube" 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,17 +67,14 @@ protected:
67 void copy_members(const ctpsfcube& app); 67 void copy_members(const ctpsfcube& app);
68 void free_members(void); 68 void free_members(void);
69 void get_parameters(void); 69 void get_parameters(void);
70 - void get_obs(void);  
71 - void set_from_cntmap(const std::string& filename);  
72 70
73 // User parameters 71 // User parameters
74 - std::string m_outfile; //!< Output PSF cube file 72 + std::string m_outcube; //!< Output PSF cube file
75 bool m_apply_edisp; //!< Apply energy dispersion? 73 bool m_apply_edisp; //!< Apply energy dispersion?
76 74
77 // Protected members 75 // Protected members
78 GObservations m_obs; //!< Observation container 76 GObservations m_obs; //!< Observation container
79 GCTAMeanPsf m_psfcube; //!< PSF cube 77 GCTAMeanPsf m_psfcube; //!< PSF cube
80 - GEbounds m_ebounds; //!< Energy boundaries  
81 }; 78 };
82 79
83 80
src/ctpsfcube/ctpsfcube.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctpsfcube.par # 7 # File: ctpsfcube.par #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Chia-Chun Lu # 10 # Author: Chia-Chun Lu #
11 # c.lu@mpi-hd.mpg.de # 11 # c.lu@mpi-hd.mpg.de #
@@ -15,11 +15,11 @@ @@ -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 exposure cube definition"
  22 +outcube, f, a, "psfcube.fits",,, "Output exposure cube file"
23 caldb, s, a, "dummy",,, "Calibration database" 23 caldb, s, a, "dummy",,, "Calibration database"
24 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 edisp, b, h, no,,, "Apply energy dispersion?" 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,6 +34,7 @@ ebinfile, f, h, &quot;NONE&quot;,,,&quot;Name of the file containing the energy bin definition&quot;
34 # 34 #
35 # Spatial binning parameters 35 # Spatial binning parameters
36 #=========================== 36 #===========================
  37 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
37 nxpix, i, a, 200,,, "Size of the X axis in pixels" 38 nxpix, i, a, 200,,, "Size of the X axis in pixels"
38 nypix, i, a, 200,,, "Size of the Y axis in pixels" 39 nypix, i, a, 200,,, "Size of the Y axis in pixels"
39 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 40 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/ctselect/ctselect.cpp
@@ -387,8 +387,7 @@ void ctselect::save(void) @@ -387,8 +387,7 @@ void ctselect::save(void)
387 void ctselect::init_members(void) 387 void ctselect::init_members(void)
388 { 388 {
389 // Initialise parameters 389 // Initialise parameters
390 - m_infile.clear();  
391 - m_outfile.clear(); 390 + m_outobs.clear();
392 m_prefix.clear(); 391 m_prefix.clear();
393 m_usepnt = false; 392 m_usepnt = false;
394 m_ra = -1.0; 393 m_ra = -1.0;
@@ -406,12 +405,6 @@ void ctselect::init_members(void) @@ -406,12 +405,6 @@ void ctselect::init_members(void)
406 m_infiles.clear(); 405 m_infiles.clear();
407 m_timemin.clear(); 406 m_timemin.clear();
408 m_timemax.clear(); 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 // Return 409 // Return
417 return; 410 return;
@@ -426,8 +419,7 @@ void ctselect::init_members(void) @@ -426,8 +419,7 @@ void ctselect::init_members(void)
426 void ctselect::copy_members(const ctselect& app) 419 void ctselect::copy_members(const ctselect& app)
427 { 420 {
428 // Copy parameters 421 // Copy parameters
429 - m_infile = app.m_infile;  
430 - m_outfile = app.m_outfile; 422 + m_outobs = app.m_outobs;
431 m_prefix = app.m_prefix; 423 m_prefix = app.m_prefix;
432 m_usepnt = app.m_usepnt; 424 m_usepnt = app.m_usepnt;
433 m_ra = app.m_ra; 425 m_ra = app.m_ra;
@@ -445,8 +437,6 @@ void ctselect::copy_members(const ctselect&amp; app) @@ -445,8 +437,6 @@ void ctselect::copy_members(const ctselect&amp; app)
445 m_infiles = app.m_infiles; 437 m_infiles = app.m_infiles;
446 m_timemin = app.m_timemin; 438 m_timemin = app.m_timemin;
447 m_timemax = app.m_timemax; 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 // Return 441 // Return
452 return; 442 return;
@@ -466,6 +456,9 @@ void ctselect::free_members(void) @@ -466,6 +456,9 @@ void ctselect::free_members(void)
466 /***********************************************************************//** 456 /***********************************************************************//**
467 * @brief Get application parameters 457 * @brief Get application parameters
468 * 458 *
  459 + * @exception GException::invalid_value
  460 + * Parameter "inobs" is required for ctselect.
  461 + *
469 * Get all task parameters from parameter file or (if required) by querying 462 * Get all task parameters from parameter file or (if required) by querying
470 * the user. Times are assumed to be in the native CTA MJD format. 463 * the user. Times are assumed to be in the native CTA MJD format.
471 * 464 *
@@ -475,40 +468,20 @@ void ctselect::free_members(void) @@ -475,40 +468,20 @@ void ctselect::free_members(void)
475 ***************************************************************************/ 468 ***************************************************************************/
476 void ctselect::get_parameters(void) 469 void ctselect::get_parameters(void)
477 { 470 {
478 - // If there are no observations in container then add a single CTA  
479 - // observation using the parameters from the parameter file 471 + // If there are no observations in container then load them via user parameters
480 if (m_obs.size() == 0) { 472 if (m_obs.size() == 0) {
481 473
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); 474 + // Throw exception if no infile is given, since this tool needs an observation
  475 + // including events
  476 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
496 477
497 - // Signal that no XML file should be used for storage  
498 - m_use_xml = false;  
499 - 478 + std::string msg = "Parameter \"inobs\" is required to be given in ctselect."
  479 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  480 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
500 } 481 }
501 -  
502 - // ... otherwise try to open as XML file  
503 - catch (GException::fits_open_error &e) {  
504 -  
505 - // Load observations from XML file  
506 - m_obs.load(m_infile);  
507 482
508 - // Signal that XML file should be used for storage  
509 - m_use_xml = true;  
510 -  
511 - } 483 + // Build observation container without response (not needed)
  484 + m_obs = get_observations(false);
512 485
513 } // endif: there was no observation in the container 486 } // endif: there was no observation in the container
514 487
@@ -537,7 +510,7 @@ void ctselect::get_parameters(void) @@ -537,7 +510,7 @@ void ctselect::get_parameters(void)
537 // Optionally read ahead parameters so that they get correctly 510 // Optionally read ahead parameters so that they get correctly
538 // dumped into the log file 511 // dumped into the log file
539 if (read_ahead()) { 512 if (read_ahead()) {
540 - m_outfile = (*this)["outfile"].filename(); 513 + m_outobs = (*this)["outobs"].filename();
541 m_prefix = (*this)["prefix"].string(); 514 m_prefix = (*this)["prefix"].string();
542 } 515 }
543 516
@@ -1048,7 +1021,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const @@ -1048,7 +1021,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const
1048 if (missing.size() > 0) { 1021 if (missing.size() > 0) {
1049 message = "The following columns are missing in the" 1022 message = "The following columns are missing in the"
1050 " \"EVENTS\" extension of input file \"" 1023 " \"EVENTS\" extension of input file \""
1051 - + m_outfile + "\": "; 1024 + + m_outobs + "\": ";
1052 for (int i = 0; i < missing.size(); ++i) { 1025 for (int i = 0; i < missing.size(); ++i) {
1053 message += "\"" + missing[i] + "\""; 1026 message += "\"" + missing[i] + "\"";
1054 if (i < missing.size()-1) { 1027 if (i < missing.size()-1) {
@@ -1060,7 +1033,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const @@ -1060,7 +1033,7 @@ std::string ctselect::check_infile(const std::string&amp; filename) const
1060 } 1033 }
1061 catch (GException::fits_hdu_not_found& e) { 1034 catch (GException::fits_hdu_not_found& e) {
1062 message = "No \"EVENTS\" extension found in input file \"" 1035 message = "No \"EVENTS\" extension found in input file \""
1063 - + m_outfile + "\"."; 1036 + + m_outobs + "\".";
1064 } 1037 }
1065 1038
1066 // Return 1039 // Return
@@ -1103,13 +1076,13 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const @@ -1103,13 +1076,13 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const
1103 void ctselect::save_fits(void) 1076 void ctselect::save_fits(void)
1104 { 1077 {
1105 // Get output filename 1078 // Get output filename
1106 - m_outfile = (*this)["outfile"].filename(); 1079 + m_outobs = (*this)["outobs"].filename();
1107 1080
1108 // Get CTA observation from observation container 1081 // Get CTA observation from observation container
1109 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]); 1082 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1110 1083
1111 // Save event list 1084 // Save event list
1112 - save_event_list(obs, m_infiles[0], m_outfile); 1085 + save_event_list(obs, m_infiles[0], m_outobs);
1113 1086
1114 // Return 1087 // Return
1115 return; 1088 return;
@@ -1129,14 +1102,14 @@ void ctselect::save_fits(void) @@ -1129,14 +1102,14 @@ void ctselect::save_fits(void)
1129 void ctselect::save_xml(void) 1102 void ctselect::save_xml(void)
1130 { 1103 {
1131 // Get output filename and prefix 1104 // Get output filename and prefix
1132 - m_outfile = (*this)["outfile"].filename(); 1105 + m_outobs = (*this)["outobs"].filename();
1133 m_prefix = (*this)["prefix"].string(); 1106 m_prefix = (*this)["prefix"].string();
1134 1107
1135 // Issue warning if output filename has no .xml suffix 1108 // Issue warning if output filename has no .xml suffix
1136 - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); 1109 + std::string suffix = gammalib::tolower(m_outobs.substr(m_outobs.length()-4,4));
1137 if (suffix != ".xml") { 1110 if (suffix != ".xml") {
1138 log << "*** WARNING: Name of observation definition output file \""+ 1111 log << "*** WARNING: Name of observation definition output file \""+
1139 - m_outfile+"\"" << std::endl; 1112 + m_outobs+"\"" << std::endl;
1140 log << "*** WARNING: does not terminate with \".xml\"." << std::endl; 1113 log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1141 log << "*** WARNING: This is not an error, but might be misleading." 1114 log << "*** WARNING: This is not an error, but might be misleading."
1142 " It is recommended" << std::endl; 1115 " It is recommended" << std::endl;
@@ -1167,7 +1140,7 @@ void ctselect::save_xml(void) @@ -1167,7 +1140,7 @@ void ctselect::save_xml(void)
1167 } // endfor: looped over observations 1140 } // endfor: looped over observations
1168 1141
1169 // Save observations in XML file 1142 // Save observations in XML file
1170 - m_obs.save(m_outfile); 1143 + m_obs.save(m_outobs);
1171 1144
1172 // Return 1145 // Return
1173 return; 1146 return;
src/ctselect/ctselect.hpp
@@ -36,7 +36,7 @@ @@ -36,7 +36,7 @@
36 36
37 /* __Definitions _________________________________________________________ */ 37 /* __Definitions _________________________________________________________ */
38 #define CTSELECT_NAME "ctselect" 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,8 +82,7 @@ protected:
82 const std::string& outfile) const; 82 const std::string& outfile) const;
83 83
84 // User parameters 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 std::string m_prefix; //!< Prefix for multiple event lists 86 std::string m_prefix; //!< Prefix for multiple event lists
88 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters 87 bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
89 double m_ra; //!< RA of ROI centre 88 double m_ra; //!< RA of ROI centre
@@ -101,8 +100,7 @@ protected: @@ -101,8 +100,7 @@ protected:
101 std::vector<std::string> m_infiles; //!< Input event filenames 100 std::vector<std::string> m_infiles; //!< Input event filenames
102 GTime m_timemin; //!< Earliest time 101 GTime m_timemin; //!< Earliest time
103 GTime m_timemax; //!< Latest time 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
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: ctselect.par # 7 # File: ctselect.par #
8 -# Version: 00-05-00 # 8 +# Version: 00-06-00 #
9 # # 9 # #
10 # Author: Jurgen Knodlseder # 10 # Author: Jurgen Knodlseder #
11 # Jurgen.Knodlseder@irap.omp.eu # 11 # Jurgen.Knodlseder@irap.omp.eu #
@@ -17,9 +17,11 @@ @@ -17,9 +17,11 @@
17 # 17 #
18 # File information 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 prefix, s, h, "selected_",,,"Prefix for event lists in observation definition file" 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 # Selection parameters 26 # Selection parameters
25 #===================== 27 #=====================
src/ctskymap/ctskymap.cpp
@@ -35,7 +35,7 @@ @@ -35,7 +35,7 @@
35 /* __ Method name definitions ____________________________________________ */ 35 /* __ Method name definitions ____________________________________________ */
36 #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)" 36 #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)"
37 #define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)" 37 #define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)"
38 - 38 +#define G_GET_PARAMETERS "ctskymap::get_parameters()"
39 /* __ Debug definitions __________________________________________________ */ 39 /* __ Debug definitions __________________________________________________ */
40 40
41 /* __ Coding definitions _________________________________________________ */ 41 /* __ Coding definitions _________________________________________________ */
@@ -262,12 +262,6 @@ void ctskymap::run(void) @@ -262,12 +262,6 @@ void ctskymap::run(void)
262 } 262 }
263 } 263 }
264 264
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 // Map events into sky map 265 // Map events into sky map
272 map_events(obs); 266 map_events(obs);
273 267
@@ -311,10 +305,10 @@ void ctskymap::save(void) @@ -311,10 +305,10 @@ void ctskymap::save(void)
311 } 305 }
312 306
313 // Get output filename 307 // Get output filename
314 - m_outfile = (*this)["outfile"].filename(); 308 + m_outmap = (*this)["outmap"].filename();
315 309
316 // Save sky map 310 // Save sky map
317 - m_skymap.save(m_outfile, clobber()); 311 + m_skymap.save(m_outmap, clobber());
318 312
319 // Return 313 // Return
320 return; 314 return;
@@ -324,6 +318,9 @@ void ctskymap::save(void) @@ -324,6 +318,9 @@ void ctskymap::save(void)
324 /***********************************************************************//** 318 /***********************************************************************//**
325 * @brief Get application parameters 319 * @brief Get application parameters
326 * 320 *
  321 + * @exception GException::invalid_value
  322 + * Parameter "inobs" is required for ctskymap.
  323 + *
327 * Get all task parameters from parameter file or (if required) by querying 324 * Get all task parameters from parameter file or (if required) by querying
328 * the user. Most parameters are only required if no observation exists so 325 * the user. Most parameters are only required if no observation exists so
329 * far in the observation container. In this case, a single CTA observation 326 * far in the observation container. In this case, a single CTA observation
@@ -332,93 +329,54 @@ void ctskymap::save(void) @@ -332,93 +329,54 @@ void ctskymap::save(void)
332 ***************************************************************************/ 329 ***************************************************************************/
333 void ctskymap::get_parameters(void) 330 void ctskymap::get_parameters(void)
334 { 331 {
335 - // If there are no observations in container then add a single CTA  
336 - // observation using the parameters from the parameter file  
337 - if (m_obs.size() == 0) { 332 + // If there are no observations in container then load them via user parameters
  333 + if (m_obs.size() == 0) {
338 334
339 - // Get name of CTA events file  
340 - m_evfile = (*this)["evfile"].filename(); 335 + // Throw exception if no infile is given, since this tool needs an observation
  336 + // including events
  337 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
341 338
342 - // Load unbinned CTA observation  
343 - GCTAObservation obs;  
344 - obs.load(m_evfile); 339 + std::string msg = "Parameter \"inobs\" is required to be given in ctskymap."
  340 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  341 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  342 + }
345 343
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(); 344 + // Build observation container without response (not needed)
  345 + m_obs = get_observations(false);
353 346
354 - } // endif: there was no observation in the container 347 + } // endif: there was no observation in the container
355 348
356 - // 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(); 349 + // Get parameters
  350 + bool usepnt = (*this)["usepnt"].boolean();
364 351
365 - // Read ahead parameters  
366 - if (read_ahead()) {  
367 - m_outfile = (*this)["outfile"].filename();  
368 - } 352 + // Check if pointing should be used
  353 + if (usepnt) {
369 354
370 - // Return  
371 - return;  
372 -} 355 + // buld map from user parameters and pointing information
  356 + m_skymap = get_map(m_obs);
373 357
  358 + } // endif: pointing used as map center
374 359
375 -/***********************************************************************//**  
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 { 360 else {
396 361
  362 + // Build cube from user parameters
  363 + m_skymap = get_map();
397 364
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 - } 365 + } // endelse: m_usepnt was false
410 366
411 - } // endelse: map centre set to pointing 367 + // Get remaining parameters
  368 + m_emin = (*this)["emin"].real();
  369 + m_emax = (*this)["emax"].real();
412 370
413 - // Create skymap  
414 - m_skymap = GSkymap(m_proj, m_coordsys, xref, yref, -m_binsz, m_binsz,  
415 - m_nxpix, m_nypix, 1); 371 + // Read ahead parameters
  372 + if (read_ahead()) {
  373 + m_outmap = (*this)["outmap"].filename();
  374 + }
416 375
417 // Return 376 // Return
418 return; 377 return;
419 } 378 }
420 379
421 -  
422 /***********************************************************************//** 380 /***********************************************************************//**
423 * @brief Map events into a sky map 381 * @brief Map events into a sky map
424 * 382 *
@@ -473,8 +431,8 @@ void ctskymap::map_events(GCTAObservation* obs) @@ -473,8 +431,8 @@ void ctskymap::map_events(GCTAObservation* obs)
473 GSkyPixel pixel = m_skymap.dir2pix(dir); 431 GSkyPixel pixel = m_skymap.dir2pix(dir);
474 432
475 // Skip if pixel is out of range 433 // 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)) { 434 + if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
  435 + pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
478 num_outside_map++; 436 num_outside_map++;
479 continue; 437 continue;
480 } 438 }
@@ -525,19 +483,10 @@ void ctskymap::map_events(GCTAObservation* obs) @@ -525,19 +483,10 @@ void ctskymap::map_events(GCTAObservation* obs)
525 void ctskymap::init_members(void) 483 void ctskymap::init_members(void)
526 { 484 {
527 // Initialise members 485 // Initialise members
528 - m_evfile.clear();  
529 - m_outfile.clear();  
530 - m_proj.clear();  
531 - m_coordsys.clear();  
532 m_obs.clear(); 486 m_obs.clear();
533 m_skymap.clear(); 487 m_skymap.clear();
534 m_emin = 0.0; 488 m_emin = 0.0;
535 m_emax = 0.0; 489 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;  
541 490
542 // Return 491 // Return
543 return; 492 return;
@@ -552,19 +501,10 @@ void ctskymap::init_members(void) @@ -552,19 +501,10 @@ void ctskymap::init_members(void)
552 void ctskymap::copy_members(const ctskymap& app) 501 void ctskymap::copy_members(const ctskymap& app)
553 { 502 {
554 // Copy attributes 503 // 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; 504 m_obs = app.m_obs;
560 m_skymap = app.m_skymap; 505 m_skymap = app.m_skymap;
561 m_emin = app.m_emin; 506 m_emin = app.m_emin;
562 m_emax = app.m_emax; 507 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;  
568 508
569 // Return 509 // Return
570 return; 510 return;
src/ctskymap/ctskymap.hpp
@@ -69,21 +69,13 @@ protected: @@ -69,21 +69,13 @@ protected:
69 void copy_members(const ctskymap& app); 69 void copy_members(const ctskymap& app);
70 void free_members(void); 70 void free_members(void);
71 void get_parameters(void); 71 void get_parameters(void);
72 - void init_map(GCTAObservation* obs);  
73 void map_events(GCTAObservation* obs); 72 void map_events(GCTAObservation* obs);
74 73
75 // User parameters 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 double m_emin; //!< Minimum energy (TeV) 76 double m_emin; //!< Minimum energy (TeV)
79 double m_emax; //!< Maximum energy (TeV) 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 // Protected members 80 // Protected members
89 GObservations m_obs; //!< Observation container 81 GObservations m_obs; //!< Observation container
src/ctskymap/ctskymap.par
@@ -17,8 +17,8 @@ @@ -17,8 +17,8 @@
17 # 17 #
18 # General parameters 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 # Event selection parameters 23 # Event selection parameters
24 #=========================== 24 #===========================
@@ -27,6 +27,7 @@ emax, r, a, 100.0,,,&quot;Maximum energy in TeV&quot; @@ -27,6 +27,7 @@ emax, r, a, 100.0,,,&quot;Maximum energy in TeV&quot;
27 # 27 #
28 # Spatial binning parameters 28 # Spatial binning parameters
29 #=========================== 29 #===========================
  30 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
30 nxpix, i, a, 200,,, "Size of the X axis in pixels" 31 nxpix, i, a, 200,,, "Size of the X axis in pixels"
31 nypix, i, a, 200,,, "Size of the Y axis in pixels" 32 nypix, i, a, 200,,, "Size of the Y axis in pixels"
32 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 33 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/cttsmap/cttsmap.cpp
@@ -233,8 +233,6 @@ void cttsmap::run(void) @@ -233,8 +233,6 @@ void cttsmap::run(void)
233 log.header1("Initialise TS map"); 233 log.header1("Initialise TS map");
234 } 234 }
235 235
236 - // Initialise maps  
237 - init_maps();  
238 236
239 // Get bins to be computed 237 // Get bins to be computed
240 // To split the computation on several jobs, 238 // To split the computation on several jobs,
@@ -369,7 +367,7 @@ void cttsmap::save(void) @@ -369,7 +367,7 @@ void cttsmap::save(void)
369 } 367 }
370 368
371 // Get output filename 369 // Get output filename
372 - m_outfile = (*this)["outfile"].filename(); 370 + m_outmap = (*this)["outmap"].filename();
373 371
374 // Create fits file instance 372 // Create fits file instance
375 GFits fitsfile = GFits(); 373 GFits fitsfile = GFits();
@@ -392,7 +390,7 @@ void cttsmap::save(void) @@ -392,7 +390,7 @@ void cttsmap::save(void)
392 } 390 }
393 391
394 // Save FITS file 392 // Save FITS file
395 - fitsfile.saveto(m_outfile,true); 393 + fitsfile.saveto(m_outmap,true);
396 394
397 // Return 395 // Return
398 return; 396 return;
@@ -411,15 +409,8 @@ void cttsmap::save(void) @@ -411,15 +409,8 @@ void cttsmap::save(void)
411 void cttsmap::init_members(void) 409 void cttsmap::init_members(void)
412 { 410 {
413 // Initialise members 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 // Initialise protected members 415 // Initialise protected members
425 m_obs.clear(); 416 m_obs.clear();
@@ -445,15 +436,8 @@ void cttsmap::init_members(void) @@ -445,15 +436,8 @@ void cttsmap::init_members(void)
445 void cttsmap::copy_members(const cttsmap& app) 436 void cttsmap::copy_members(const cttsmap& app)
446 { 437 {
447 // Copy attributes 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 // Copy protected members 442 // Copy protected members
459 m_binmin = app.m_binmin; 443 m_binmin = app.m_binmin;
@@ -486,6 +470,9 @@ void cttsmap::free_members(void) @@ -486,6 +470,9 @@ void cttsmap::free_members(void)
486 /***********************************************************************//** 470 /***********************************************************************//**
487 * @brief Get application parameters 471 * @brief Get application parameters
488 * 472 *
  473 + * @exception GException::invalid_value
  474 + * Parameter "inobs" is required for cttsmap.
  475 + *
489 * @exception GException::invalid_value 476 * @exception GException::invalid_value
490 * Test source not found or no RA/DEC parameters found for test 477 * Test source not found or no RA/DEC parameters found for test
491 * source. 478 * source.
@@ -498,42 +485,20 @@ void cttsmap::free_members(void) @@ -498,42 +485,20 @@ void cttsmap::free_members(void)
498 ***************************************************************************/ 485 ***************************************************************************/
499 void cttsmap::get_parameters(void) 486 void cttsmap::get_parameters(void)
500 { 487 {
501 - // If there are no observations in container then add a single CTA  
502 - // observation using the parameters from the parameter file 488 + // If there are no observations in container then load them via user parameters
503 if (m_obs.size() == 0) { 489 if (m_obs.size() == 0) {
504 490
505 - // Allocate CTA observation  
506 - GCTAObservation obs;  
507 -  
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); 491 + // Throw exception if no infile is given, since this tool needs an observation
  492 + // including events
  493 + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") {
522 494
  495 + std::string msg = "Parameter \"inobs\" is required to be given in cttsmap."
  496 + "Specify a vaild observation definition (XML or FITS) file to proceed";
  497 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
523 } 498 }
524 499
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 500 + // Build observation container
  501 + m_obs = get_observations();
537 502
538 } // endif: there was no observation in the container 503 } // endif: there was no observation in the container
539 504
@@ -542,7 +507,7 @@ void cttsmap::get_parameters(void) @@ -542,7 +507,7 @@ void cttsmap::get_parameters(void)
542 if (m_obs.models().size() == 0) { 507 if (m_obs.models().size() == 0) {
543 508
544 // Get models XML filename 509 // Get models XML filename
545 - std::string filename = (*this)["srcmdl"].filename(); 510 + std::string filename = (*this)["inmodel"].filename();
546 511
547 // Setup models for optimizing. 512 // Setup models for optimizing.
548 m_obs.models(GModels(filename)); 513 m_obs.models(GModels(filename));
@@ -579,14 +544,27 @@ void cttsmap::get_parameters(void) @@ -579,14 +544,27 @@ void cttsmap::get_parameters(void)
579 (*m_testsource)["DEC"].fix(); 544 (*m_testsource)["DEC"].fix();
580 } 545 }
581 546
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(); 547 +
  548 + // Get sky map binning from user parameters
  549 + // and check if pointing should be used
  550 + GSkymap map;
  551 + bool usepnt = (*this)["usepnt"].boolean();
  552 + if (usepnt) {
  553 +
  554 + // build cube from user parameters and pointing information
  555 + map = get_map(m_obs);
  556 +
  557 + } // endif: pointing used as map center
  558 +
  559 + else {
  560 +
  561 + // Build cube from user parameters
  562 + map = get_map();
  563 +
  564 + } // endelse: m_usepnt was false
  565 +
  566 + // Initialise maps from user parameters
  567 + init_maps(map);
590 568
591 // Get optional splitting parameters 569 // Get optional splitting parameters
592 m_binmin = (*this)["binmin"].integer(); 570 m_binmin = (*this)["binmin"].integer();
@@ -596,7 +574,7 @@ void cttsmap::get_parameters(void) @@ -596,7 +574,7 @@ void cttsmap::get_parameters(void)
596 // Optionally read ahead parameters so that they get correctly 574 // Optionally read ahead parameters so that they get correctly
597 // dumped into the log file 575 // dumped into the log file
598 if (read_ahead()) { 576 if (read_ahead()) {
599 - m_outfile = (*this)["outfile"].filename(); 577 + m_outmap = (*this)["outmap"].filename();
600 } 578 }
601 579
602 // Return 580 // Return
@@ -609,23 +587,19 @@ void cttsmap::get_parameters(void) @@ -609,23 +587,19 @@ void cttsmap::get_parameters(void)
609 * 587 *
610 * Initialises skymaps that will contain map information. 588 * Initialises skymaps that will contain map information.
611 ***************************************************************************/ 589 ***************************************************************************/
612 -void cttsmap::init_maps(void) 590 +void cttsmap::init_maps(const GSkymap& map)
613 { 591 {
614 // Initialise map information 592 // Initialise map information
615 m_tsmap.clear(); 593 m_tsmap.clear();
616 594
617 // Create skymap 595 // 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); 596 + m_tsmap = GSkymap(map);
621 597
622 // Initialise map information 598 // Initialise map information
623 m_statusmap.clear(); 599 m_statusmap.clear();
624 600
625 // Create status map 601 // 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); 602 + m_statusmap = GSkymap(map);
629 603
630 // Initialise maps of free parameters 604 // Initialise maps of free parameters
631 if (m_testsource != NULL) { 605 if (m_testsource != NULL) {
@@ -640,7 +614,7 @@ void cttsmap::init_maps(void) @@ -640,7 +614,7 @@ void cttsmap::init_maps(void)
640 614
641 // Add sky map for free parameters. Note that push_back will 615 // Add sky map for free parameters. Note that push_back will
642 // create a copy. 616 // create a copy.
643 - m_maps.push_back(m_tsmap); 617 + m_maps.push_back(map);
644 618
645 // Store parameter name 619 // Store parameter name
646 m_mapnames.push_back((*m_testsource)[i].name()); 620 m_mapnames.push_back((*m_testsource)[i].name());
src/cttsmap/cttsmap.hpp
@@ -34,7 +34,7 @@ @@ -34,7 +34,7 @@
34 34
35 /* __Definitions _________________________________________________________ */ 35 /* __Definitions _________________________________________________________ */
36 #define CTTSMAP_NAME "cttsmap" 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,19 +81,11 @@ protected: @@ -81,19 +81,11 @@ protected:
81 void copy_members(const cttsmap& app); 81 void copy_members(const cttsmap& app);
82 void free_members(void); 82 void free_members(void);
83 void get_parameters(void); 83 void get_parameters(void);
84 - void init_maps(void); 84 + void init_maps(const GSkymap& map);
85 85
86 // User parameters 86 // User parameters
87 - std::string m_infile; //!< Input file  
88 std::string m_srcname; //!< Name of source which is moved around 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 // Parameters to control speed and job splitting 90 // Parameters to control speed and job splitting
99 int m_binmin; //!< Map bin number from which computation should start 91 int m_binmin; //!< Map bin number from which computation should start
src/cttsmap/cttsmap.par
@@ -5,7 +5,7 @@ @@ -5,7 +5,7 @@
5 # ------------------------------------------------------------------------- # 5 # ------------------------------------------------------------------------- #
6 # # 6 # #
7 # File: cttsmap.par # 7 # File: cttsmap.par #
8 -# Version: 00-01-00 # 8 +# Version: 00-02-00 #
9 # # 9 # #
10 # Author: Michael Mayer # 10 # Author: Michael Mayer #
11 # michael.mayer@desy.de # 11 # michael.mayer@desy.de #
@@ -17,16 +17,17 @@ @@ -17,16 +17,17 @@
17 # 17 #
18 # General parameters 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 srcname, s, a, "TestSource",,,"Test source" 22 srcname, s, a, "TestSource",,,"Test source"
23 caldb, s, a, "dummy",,, "Calibration database" 23 caldb, s, a, "dummy",,, "Calibration database"
24 irf, s, a, "cta_dummy_irf",,, "Instrument response function" 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 -outfile, f, a, "tsmap.fits",,, "Output counts cube" 25 +outmap, f, a, "tsmap.fits",,, "Output counts cube"
26 26
27 # 27 #
28 # Spatial binning parameters 28 # Spatial binning parameters
29 #=========================== 29 #===========================
  30 +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
30 nxpix, i, a, 200,,, "Size of the X axis in pixels" 31 nxpix, i, a, 200,,, "Size of the X axis in pixels"
31 nypix, i, a, 200,,, "Size of the Y axis in pixels" 32 nypix, i, a, 200,,, "Size of the Y axis in pixels"
32 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" 33 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
src/support/ctool.cpp
@@ -48,6 +48,8 @@ @@ -48,6 +48,8 @@
48 48
49 /* __ Method name definitions ____________________________________________ */ 49 /* __ Method name definitions ____________________________________________ */
50 #define G_GET_EBOUNDS "ctool::get_ebounds()" 50 #define G_GET_EBOUNDS "ctool::get_ebounds()"
  51 +#define G_SET_FROM_CNTMAP "ctool::set_from_cntmap(std::string&)"
  52 +#define G_GET_CUBE "ctool::get_cube(GObservations&)"
51 53
52 /* __ Debug definitions __________________________________________________ */ 54 /* __ Debug definitions __________________________________________________ */
53 55
@@ -226,6 +228,12 @@ void ctool::init_members(void) @@ -226,6 +228,12 @@ void ctool::init_members(void)
226 { 228 {
227 // Initialise members 229 // Initialise members
228 m_read_ahead = false; 230 m_read_ahead = false;
  231 + m_use_xml = false;
  232 +
  233 + // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
  234 + // which is defined in GCTALib.hpp. This is somehow a kluge. We need
  235 + // a better mechanism to implement the CTA reference MJD.
  236 + m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
229 237
230 // Set logger properties 238 // Set logger properties
231 log.date(true); 239 log.date(true);
@@ -244,6 +252,8 @@ void ctool::copy_members(const ctool&amp; app) @@ -244,6 +252,8 @@ void ctool::copy_members(const ctool&amp; app)
244 { 252 {
245 // Copy members 253 // Copy members
246 m_read_ahead = app.m_read_ahead; 254 m_read_ahead = app.m_read_ahead;
  255 + m_use_xml = app.m_use_xml;
  256 + m_cta_ref = app.m_cta_ref;
247 257
248 // Return 258 // Return
249 return; 259 return;
@@ -342,6 +352,412 @@ GEbounds ctool::get_ebounds(void) @@ -342,6 +352,412 @@ GEbounds ctool::get_ebounds(void)
342 } 352 }
343 353
344 354
  355 +
  356 +
  357 +/***********************************************************************//**
  358 + * @brief Get a GCTAEventCube from file
  359 + *
  360 + * @param[in] filename FITS file of count cube.
  361 + * @exception GException::invalid_value
  362 + * Counts map not of WCS type
  363 + *
  364 + * Gets an cta event cube from file. This function is used in many ctools
  365 + * to obtain a skymap definition from file
  366 + *
  367 + ***************************************************************************/
  368 +GCTAEventCube ctool::set_from_cntmap(const std::string filename)
  369 +{
  370 + // Allocate CTA observation
  371 + GCTAEventCube cube;
  372 +
  373 + // Load counts map in CTA observation
  374 + cube.load(filename);
  375 +
  376 + // Check for valid WCS coordinate system
  377 + const GWcs* wcs = dynamic_cast<const GWcs*>(cube.map().projection());
  378 + if (wcs == NULL) {
  379 +
  380 + std::string msg = "Counts map project is not of WCS type.";
  381 + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  382 +
  383 + } // endif: WCS projection was not valid
  384 +
  385 + // Return cube
  386 + return cube;
  387 +}
  388 +
  389 +
  390 +
  391 +/***********************************************************************//**
  392 + * @brief Get observation container
  393 + *
  394 + * @param[in] get_response flag whether response information should be loaded
  395 + *
  396 + * Get an observation container according to the user parameters. The method
  397 + * supports loading of a individual FITS file, an observation definition file in XML
  398 + * format. If the input filename is empty, parameters are read to build a CTA observation
  399 + * from scratch.
  400 + *
  401 + ***************************************************************************/
  402 +GObservations ctool::get_observations(const bool& get_response)
  403 +{
  404 + // Initialise empty observation container
  405 + GObservations obs;
  406 +
  407 + // Get the filename from the input parameters
  408 + std::string filename = (*this)["inobs"].filename();
  409 +
  410 + // If no observation definition file has been specified then read all
  411 + // parameters that are necessary to create an observation from scratch
  412 + // (see method setup_obs)
  413 + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) {
  414 +
  415 + // Setup a new CTA observation
  416 + GCTAObservation cta_obs = setup_obs();
  417 +
  418 + // Get response if required
  419 + if (get_response) {
  420 +
  421 + // Set response
  422 + set_obs_response(&cta_obs);
  423 +
  424 + } // endif: response was required
  425 +
  426 + // Append observation to container
  427 + obs.append(cta_obs);
  428 +
  429 + } // endif: filename was "NONE" or ""
  430 +
  431 + else {
  432 +
  433 + // Try first to open as FITS file
  434 + try {
  435 +
  436 + // Allocate empyt CTA observation
  437 + GCTAObservation cta_obs;
  438 +
  439 + // Load data
  440 + cta_obs.load(filename);
  441 +
  442 + // Get response if required
  443 + if (get_response) {
  444 +
  445 + // Set response
  446 + set_obs_response(&cta_obs);
  447 +
  448 + } // endif: response was required
  449 +
  450 + // Append observation to container
  451 + obs.append(cta_obs);
  452 +
  453 + // Signal that no XML file should be used for storage
  454 + m_use_xml = false;
  455 +
  456 + }
  457 + // ... otherwise try to open as XML file
  458 + catch (GException::fits_open_error &e) {
  459 +
  460 + // Load observations from XML file
  461 + obs.load(filename);
  462 +
  463 + // Get response if required
  464 + if (get_response) {
  465 +
  466 + // For all observations that have no response, set the response
  467 + // from the task parameters
  468 + set_response(obs);
  469 +
  470 + } // endif: response was required
  471 +
  472 + // Signal that XML file should be used for storage
  473 + m_use_xml = true;
  474 +
  475 + } // endcatch: file was an XML file
  476 + }
  477 +
  478 + // Return observation container
  479 + return obs;
  480 +
  481 +}
  482 +
  483 +/***********************************************************************//**
  484 + * @brief Build a CTA observation with an empty event list
  485 + *
  486 + * Creates an empty CTA observation from user parameters. An empty event list
  487 + * including RoI, GTi and energy boundary information is attached to the observation.
  488 + * The method also sets the pointing direction using the ra and dec parameter, the ROI based on
  489 + * ra, dec and rad, a single GTI based on tmin and tmax, and a
  490 + * single energy boundary based on emin and emax. The method furthermore
  491 + * sets the ontime, livetime and deadtime correction factor.
  492 + *
  493 + ***************************************************************************/
  494 +GCTAObservation ctool::setup_obs(void)
  495 +{
  496 + //Get CTA observation parameters
  497 + double ra = (*this)["ra"].real();
  498 + double dec = (*this)["dec"].real();
  499 + double rad = (*this)["rad"].real();
  500 + double deadc = (*this)["deadc"].real();
  501 + double t_min = (*this)["tmin"].real();
  502 + double t_max = (*this)["tmax"].real();
  503 + double e_min = (*this)["emin"].real();
  504 + double e_max = (*this)["emax"].real();
  505 +
  506 + // Allocate CTA observation and empty event list
  507 + GCTAObservation obs;
  508 + GCTAEventList list;
  509 +
  510 + // Set pointing direction
  511 + GCTAPointing pnt;
  512 + GSkyDir skydir;
  513 + skydir.radec_deg(ra, dec);
  514 + pnt.dir(skydir);
  515 +
  516 + // Set ROI
  517 + GCTAInstDir instdir(skydir);
  518 + GCTARoi roi(instdir, rad);
  519 +
  520 + // Set GTI
  521 + GGti gti(m_cta_ref);
  522 + GTime tstart;
  523 + GTime tstop;
  524 + tstart.set(t_min, m_cta_ref);
  525 + tstop.set(t_max, m_cta_ref);
  526 + gti.append(tstart, tstop);
  527 +
  528 + // Set energy boundaries
  529 + GEbounds ebounds;
  530 + GEnergy emin;
  531 + GEnergy emax;
  532 + emin.TeV(e_min);
  533 + emax.TeV(e_max);
  534 + ebounds.append(emin, emax);
  535 +
  536 + // Set CTA event list attributes
  537 + list.roi(roi);
  538 + list.gti(gti);
  539 + list.ebounds(ebounds);
  540 +
  541 + // Attach event list to CTA observation
  542 + obs.events(list);
  543 +
  544 + // Set observation ontime, livetime and deadtime correction factor
  545 + obs.pointing(pnt);
  546 + obs.ontime(gti.ontime());
  547 + obs.livetime(gti.ontime()*deadc);
  548 + obs.deadc(deadc);
  549 +
  550 + // Return CTA observation
  551 + return obs;
  552 +
  553 +}
  554 +
  555 +GSkyDir ctool::get_pointing(const GObservations& obs)
  556 +{
  557 + // Initialise values
  558 + double xref = 0.0;
  559 + double yref = 0.0;
  560 + int n_pnt = 0;
  561 + bool found = false;
  562 +
  563 + // Loop over all observations to compute mean pointings
  564 + for (int i = 0; i < obs.size(); ++i) {
  565 + const GCTAObservation* cta_obs = dynamic_cast<const GCTAObservation*>(obs[i]);
  566 + if (cta_obs != NULL) {
  567 +
  568 + // get pointing
  569 + const GCTAPointing& pnt = cta_obs->pointing();
  570 +
  571 + // Add to coordinates
  572 + xref += pnt.dir().ra_deg();
  573 + yref += pnt.dir().dec_deg();
  574 + n_pnt += 1;
  575 +
  576 + // Signal that pointing has been found
  577 + found = true;
  578 +
  579 + } // endif: cta observation was valid
  580 +
  581 + } // endfor: loop over alle observations
  582 +
  583 + // Signal if no pointing is found
  584 + if (!found) {
  585 + std::string msg = "No valid CTA observation has been found in "
  586 + "observation list, hence no pointing information "
  587 + "could be extracted. Use the \"usepnt=no\" "
  588 + "option and specify pointing explicitly";
  589 + throw GException::invalid_value(G_GET_CUBE, msg);
  590 + }
  591 +
  592 + // Calculate mean coordinates
  593 + double ra = xref / double (n_pnt);
  594 + double dec = yref / double (n_pnt);
  595 +
  596 + // Initialise sky dir
  597 + GSkyDir dir = GSkyDir();
  598 + dir.radec_deg(ra,dec);
  599 +
  600 + // Return sky direction
  601 + return dir;
  602 +}
  603 +
  604 +/***********************************************************************//**
  605 + * @brief Build a GCTAEventCube based on user parameters
  606 + *
  607 + * Creates an empty CTA event cube from input user parameters. It appends a dummy
  608 + * GTi to the event cube.
  609 + *
  610 + ***************************************************************************/
  611 +GCTAEventCube ctool::get_cube(void)
  612 +{
  613 + // Get cube parameters
  614 + double xref = (*this)["xref"].real();
  615 + double yref = (*this)["yref"].real();
  616 + std::string proj = (*this)["proj"].string();
  617 + std::string coordsys = (*this)["coordsys"].string();
  618 + double binsz = (*this)["binsz"].real();
  619 + int nxpix = (*this)["nxpix"].integer();
  620 + int nypix = (*this)["nypix"].integer();
  621 +
  622 + // Set energy boundaries
  623 + GEbounds ebounds = get_ebounds();
  624 +
  625 + // Set dummy GTI that is needed for event cube creation
  626 + GGti gti;
  627 + gti.append(GTime(0.0), GTime(0.1234));
  628 +
  629 + // Initialise sky map
  630 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, ebounds.size());
  631 +
  632 + // Build event cube from given information
  633 + GCTAEventCube cube = GCTAEventCube(map,ebounds,gti);
  634 +
  635 + // Return cube
  636 + return cube;
  637 +
  638 +}
  639 +
  640 +/***********************************************************************//**
  641 + * @brief Build a GCTAEventCube based on user parameters
  642 + *
  643 + * @param[in] obs observation container
  644 + *
  645 + * Creates an empty CTA event cube from input user parameters. It appends a dummy
  646 + * GTi to the event cube. The input parameter is used to extract the pointing information
  647 + * from the given observations. This function is used if a ctool has the input parameter
  648 + * usepnt=true.
  649 + *
  650 + ***************************************************************************/
  651 +GCTAEventCube ctool::get_cube(const GObservations& obs)
  652 +{
  653 + // Read ahead map information
  654 + std::string proj = (*this)["proj"].string();
  655 + std::string coordsys = (*this)["coordsys"].string();
  656 +
  657 + // Initialise reference coordinates
  658 + double xref = 0.0;
  659 + double yref = 0.0;
  660 +
  661 + // Get pointing from observations
  662 + GSkyDir pnt = get_pointing(obs);
  663 +
  664 + // Use correct coordinate system
  665 + if (gammalib::toupper(coordsys) == "GAL") {
  666 + xref = pnt.l_deg();
  667 + yref = pnt.b_deg();
  668 + }
  669 + else {
  670 + xref = pnt.ra_deg();
  671 + yref = pnt.dec_deg();
  672 + }
  673 +
  674 + // Read further binning parameters
  675 + double binsz = (*this)["binsz"].real();
  676 + int nxpix = (*this)["nxpix"].integer();
  677 + int nypix = (*this)["nypix"].integer();
  678 +
  679 + // Set energy boundaries
  680 + GEbounds ebounds = get_ebounds();
  681 +
  682 + // Set dummy GTI that is needed for event cube creation
  683 + GGti gti;
  684 + gti.append(GTime(0.0), GTime(0.1234));
  685 +
  686 + // Build skymap from user parameters
  687 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, ebounds.size());
  688 +
  689 + // Initialise event cube from all parameters
  690 + GCTAEventCube cube = GCTAEventCube(map,ebounds,gti);
  691 +
  692 + // Return event cube
  693 + return cube;
  694 +
  695 +}
  696 +
  697 +/***********************************************************************//**
  698 + * @brief Build a skymap based on user parameters
  699 + *
  700 + * Creates an empty skymap from input user parameters.
  701 + *
  702 + ***************************************************************************/
  703 +GSkymap ctool::get_map()
  704 +{
  705 + // Read map information
  706 + double xref = (*this)["xref"].real();
  707 + double yref = (*this)["yref"].real();
  708 + std::string proj = (*this)["proj"].string();
  709 + std::string coordsys = (*this)["coordsys"].string();
  710 + double binsz = (*this)["binsz"].real();
  711 + int nxpix = (*this)["nxpix"].integer();
  712 + int nypix = (*this)["nypix"].integer();
  713 +
  714 + // Initialise sky map
  715 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, 1);
  716 +
  717 + // Return sky map
  718 + return map;
  719 +}
  720 +
  721 +/***********************************************************************//**
  722 + * @brief Build a skymap based on user parameters
  723 + *
  724 + * Creates an empty skymap from input user parameters.
  725 + *
  726 + ***************************************************************************/
  727 +GSkymap ctool::get_map(const GObservations& obs)
  728 +{
  729 + // Read ahead map information
  730 + std::string proj = (*this)["proj"].string();
  731 + std::string coordsys = (*this)["coordsys"].string();
  732 +
  733 + // Initialise reference coordinates
  734 + double xref = 0.0;
  735 + double yref = 0.0;
  736 +
  737 + // Get pointing from observations
  738 + GSkyDir pnt = get_pointing(obs);
  739 +
  740 + // Use correct coordinate system
  741 + if (gammalib::toupper(coordsys) == "GAL") {
  742 + xref = pnt.l_deg();
  743 + yref = pnt.b_deg();
  744 + }
  745 + else {
  746 + xref = pnt.ra_deg();
  747 + yref = pnt.dec_deg();
  748 + }
  749 +
  750 + double binsz = (*this)["binsz"].real();
  751 + int nxpix = (*this)["nxpix"].integer();
  752 + int nypix = (*this)["nypix"].integer();
  753 +
  754 + // Initialise sky map
  755 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, 1);
  756 +
  757 + // Return sky map
  758 + return map;
  759 +}
  760 +
345 /***********************************************************************//** 761 /***********************************************************************//**
346 * @brief Set response for all CTA observations in container 762 * @brief Set response for all CTA observations in container
347 * 763 *
src/support/ctool.hpp
@@ -67,12 +67,23 @@ protected: @@ -67,12 +67,23 @@ protected:
67 void free_members(void); 67 void free_members(void);
68 const bool& read_ahead(void) const; 68 const bool& read_ahead(void) const;
69 GEbounds get_ebounds(void); 69 GEbounds get_ebounds(void);
  70 + GObservations get_observations(const bool& get_response=true);
  71 + GCTAObservation setup_obs(void);
  72 + GSkymap get_map(void);
  73 + GCTAEventCube get_cube(void);
  74 + GSkymap get_map(const GObservations& obs);
  75 + GCTAEventCube get_cube(const GObservations& obs);
  76 + GCTAEventCube set_from_cntmap(const std::string filename);
  77 + GSkyDir get_pointing(const GObservations& obs);
  78 +
70 void set_response(GObservations& obs); 79 void set_response(GObservations& obs);
71 void set_obs_response(GCTAObservation* obs); 80 void set_obs_response(GCTAObservation* obs);
72 size_t get_current_rss(void); 81 size_t get_current_rss(void);
73 82
74 // Protected members 83 // Protected members
75 bool m_read_ahead; //!< Read ahead parameters 84 bool m_read_ahead; //!< Read ahead parameters
  85 + bool m_use_xml; //!< Use XML file instead of FITS file for observations
  86 + GTimeReference m_cta_ref; //!< CTA time reference
76 }; 87 };
77 88
78 89
test/test_cscripts.sh.in
@@ -52,7 +52,7 @@ export PFILES=pfiles @@ -52,7 +52,7 @@ export PFILES=pfiles
52 $ECHO @ECHO_N@ "Test cspull: " 52 $ECHO @ECHO_N@ "Test cspull: "
53 chmod +x @abs_top_srcdir@/scripts/cspull.py 53 chmod +x @abs_top_srcdir@/scripts/cspull.py
54 @abs_top_srcdir@/scripts/cspull.py \ 54 @abs_top_srcdir@/scripts/cspull.py \
55 - srcmdl="data/crab.xml" \ 55 + inmodel="data/crab.xml" \
56 outfile="pull.dat" \ 56 outfile="pull.dat" \
57 ntrials="3" \ 57 ntrials="3" \
58 caldb="irf" \ 58 caldb="irf" \
@@ -138,9 +138,9 @@ $ECHO &quot; ok&quot; @@ -138,9 +138,9 @@ $ECHO &quot; ok&quot;
138 $ECHO @ECHO_N@ "Test csresmap: " 138 $ECHO @ECHO_N@ "Test csresmap: "
139 chmod +x @abs_top_srcdir@/scripts/csresmap.py 139 chmod +x @abs_top_srcdir@/scripts/csresmap.py
140 @abs_top_srcdir@/scripts/csresmap.py \ 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 caldb="irf" \ 144 caldb="irf" \
145 irf="cta_dummy_irf" \ 145 irf="cta_dummy_irf" \
146 nxpix=50 \ 146 nxpix=50 \
test/test_ctbin.py
@@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 """ 64 """
65 # Set-up ctbin 65 # Set-up ctbin
66 bin = ctools.ctbin() 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 bin["ebinalg"].string("LOG") 69 bin["ebinalg"].string("LOG")
70 bin["emin"].real(0.1) 70 bin["emin"].real(0.1)
71 bin["emax"].real(100.0) 71 bin["emax"].real(100.0)
test/test_ctbkgcube.py
@@ -40,6 +40,8 @@ class Test(gammalib.GPythonTestSuite): @@ -40,6 +40,8 @@ class Test(gammalib.GPythonTestSuite):
40 # Set members 40 # Set members
41 self.events_name = "data/crab_events.fits" 41 self.events_name = "data/crab_events.fits"
42 self.bkg_model = "data/crab.xml" 42 self.bkg_model = "data/crab.xml"
  43 + self.caldb = "irf"
  44 + self.irf = "cta_dummy_irf"
43 45
44 # Return 46 # Return
45 return 47 return
@@ -65,11 +67,13 @@ class Test(gammalib.GPythonTestSuite): @@ -65,11 +67,13 @@ class Test(gammalib.GPythonTestSuite):
65 """ 67 """
66 # Set-up ctbkgcube 68 # Set-up ctbkgcube
67 bkgcube = ctools.ctbkgcube() 69 bkgcube = ctools.ctbkgcube()
68 - bkgcube["infile"].filename(self.events_name) 70 + bkgcube["inobs"].filename(self.events_name)
69 bkgcube["inmodel"].filename(self.bkg_model) 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 bkgcube["outmodel"].filename("bkgcube.xml") 74 bkgcube["outmodel"].filename("bkgcube.xml")
  75 + bkgcube["caldb"].string(self.caldb)
  76 + bkgcube["irf"].string(self.irf)
73 bkgcube["ebinalg"].string("LOG") 77 bkgcube["ebinalg"].string("LOG")
74 bkgcube["emin"].real(0.1) 78 bkgcube["emin"].real(0.1)
75 bkgcube["emax"].real(100.0) 79 bkgcube["emax"].real(100.0)
test/test_ctbutterfly.py
@@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite): @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite):
67 """ 67 """
68 # Set-up ctbutterfly 68 # Set-up ctbutterfly
69 butterfly = ctools.ctbutterfly() 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 butterfly["srcname"].string("Crab") 72 butterfly["srcname"].string("Crab")
73 butterfly["caldb"].string(self.caldb) 73 butterfly["caldb"].string(self.caldb)
74 butterfly["irf"].string(self.irf) 74 butterfly["irf"].string(self.irf)
test/test_ctcubemask.py
@@ -65,9 +65,9 @@ class Test(gammalib.GPythonTestSuite): @@ -65,9 +65,9 @@ class Test(gammalib.GPythonTestSuite):
65 """ 65 """
66 # Set-up ctcubemask 66 # Set-up ctcubemask
67 mask = ctools.ctcubemask() 67 mask = ctools.ctcubemask()
68 - mask["infile"].filename(self.cntmap_name) 68 + mask["inobs"].filename(self.cntmap_name)
69 mask["regfile"].filename(self.regfile) 69 mask["regfile"].filename(self.regfile)
70 - mask["outfile"].filename("filtered_cntmap.fits") 70 + mask["outcube"].filename("filtered_cntmap.fits")
71 mask["ra"].real(83.63) 71 mask["ra"].real(83.63)
72 mask["dec"].real(22.01) 72 mask["dec"].real(22.01)
73 mask["rad"].real(2.0) 73 mask["rad"].real(2.0)
test/test_ctexpcube.py
@@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite): @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite):
66 """ 66 """
67 # Set-up ctexpcube 67 # Set-up ctexpcube
68 expcube = ctools.ctexpcube() 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 expcube["caldb"].string(self.caldb) 72 expcube["caldb"].string(self.caldb)
73 expcube["irf"].string(self.irf) 73 expcube["irf"].string(self.irf)
74 expcube["ebinalg"].string("LOG") 74 expcube["ebinalg"].string("LOG")
test/test_ctlike.py
@@ -67,11 +67,11 @@ class Test(gammalib.GPythonTestSuite): @@ -67,11 +67,11 @@ class Test(gammalib.GPythonTestSuite):
67 """ 67 """
68 # Set-up ctlike 68 # Set-up ctlike
69 like = ctools.ctlike() 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 like["caldb"].string(self.caldb) 72 like["caldb"].string(self.caldb)
73 like["irf"].string(self.irf) 73 like["irf"].string(self.irf)
74 - like["outmdl"].filename("result.xml") 74 + like["outmodel"].filename("result.xml")
75 75
76 # Run tool 76 # Run tool
77 self.test_try("Run ctlike") 77 self.test_try("Run ctlike")
test/test_ctmodel.py
@@ -66,12 +66,13 @@ class Test(gammalib.GPythonTestSuite): @@ -66,12 +66,13 @@ class Test(gammalib.GPythonTestSuite):
66 """ 66 """
67 # Set-up ctmodel from scratch 67 # Set-up ctmodel from scratch
68 model = ctools.ctmodel() 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 model["caldb"].string(self.caldb) 73 model["caldb"].string(self.caldb)
74 model["irf"].string(self.irf) 74 model["irf"].string(self.irf)
  75 + model["rad"].real(5.0)
75 model["ra"].real(83.63) 76 model["ra"].real(83.63)
76 model["dec"].real(22.01) 77 model["dec"].real(22.01)
77 model["tmin"].real(0.0) 78 model["tmin"].real(0.0)
test/test_ctobssim.py
@@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite): @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite):
67 """ 67 """
68 # Set-up ctobssim 68 # Set-up ctobssim
69 sim = ctools.ctobssim() 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 sim["caldb"].string(self.caldb) 72 sim["caldb"].string(self.caldb)
73 sim["irf"].string(self.irf) 73 sim["irf"].string(self.irf)
74 sim["ra"].real(83.63) 74 sim["ra"].real(83.63)
@@ -117,8 +117,8 @@ class Test(gammalib.GPythonTestSuite): @@ -117,8 +117,8 @@ class Test(gammalib.GPythonTestSuite):
117 117
118 # Set-up ctobssim 118 # Set-up ctobssim
119 sim = ctools.ctobssim(obs) 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 # Run tool 122 # Run tool
123 self.test_try("Run ctobssim") 123 self.test_try("Run ctobssim")
124 try: 124 try:
test/test_ctools.sh
@@ -67,10 +67,11 @@ export PFILES=pfiles @@ -67,10 +67,11 @@ export PFILES=pfiles
67 # Test ctobssim 67 # Test ctobssim
68 # ============= 68 # =============
69 $ECHO -n "Test ctobssim: " 69 $ECHO -n "Test ctobssim: "
  70 +
70 # 71 #
71 # Run 1 72 # Run 1
72 -$ctobssim infile="data/crab.xml" \  
73 - outfile="events.fits" \ 73 +$ctobssim inmodel="data/crab.xml" \
  74 + outevents="events.fits" \
74 caldb="irf" \ 75 caldb="irf" \
75 irf="cta_dummy_irf" \ 76 irf="cta_dummy_irf" \
76 ra=83.63 \ 77 ra=83.63 \
@@ -90,8 +91,8 @@ else @@ -90,8 +91,8 @@ else
90 fi 91 fi
91 # 92 #
92 # Run 2 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 caldb="irf" \ 96 caldb="irf" \
96 irf="cta_dummy_irf" \ 97 irf="cta_dummy_irf" \
97 ra=84.17263 \ 98 ra=84.17263 \
@@ -116,8 +117,8 @@ $ECHO &quot; ok&quot; @@ -116,8 +117,8 @@ $ECHO &quot; ok&quot;
116 # Test ctskymap 117 # Test ctskymap
117 # ============= 118 # =============
118 $ECHO -n "Test ctskymap: " 119 $ECHO -n "Test ctskymap: "
119 -$ctskymap evfile="events.fits" \  
120 - outfile="skymap.fits" \ 120 +$ctskymap inobs="events.fits" \
  121 + outmap="skymap.fits" \
121 emin=0.1 \ 122 emin=0.1 \
122 emax=100.0 \ 123 emax=100.0 \
123 nxpix=200 \ 124 nxpix=200 \
@@ -144,8 +145,8 @@ $ECHO &quot; ok&quot; @@ -144,8 +145,8 @@ $ECHO &quot; ok&quot;
144 $ECHO -n "Test ctbin: " 145 $ECHO -n "Test ctbin: "
145 # 146 #
146 # Run 1 147 # Run 1
147 -$ctbin evfile="events.fits" \  
148 - outfile="cntmap.fits" \ 148 +$ctbin inobs="events.fits" \
  149 + outcube="cntmap.fits" \
149 emin=0.1 \ 150 emin=0.1 \
150 emax=100.0 \ 151 emax=100.0 \
151 enumbins=20 \ 152 enumbins=20 \
@@ -167,8 +168,8 @@ else @@ -167,8 +168,8 @@ else
167 fi 168 fi
168 # 169 #
169 # Run 2 170 # Run 2
170 -$ctbin evfile="events2.fits" \  
171 - outfile="cntmap2.fits" \ 171 +$ctbin inobs="events2.fits" \
  172 + outcube="cntmap2.fits" \
172 emin=0.3 \ 173 emin=0.3 \
173 emax=8.0 \ 174 emax=8.0 \
174 enumbins=10 \ 175 enumbins=10 \
@@ -191,18 +192,19 @@ fi @@ -191,18 +192,19 @@ fi
191 $ECHO " ok" 192 $ECHO " ok"
192 193
193 194
  195 +
194 # 196 #
195 # Test ctmodel 197 # Test ctmodel
196 # ============ 198 # ============
197 $ECHO -n "Test ctmodel: " 199 $ECHO -n "Test ctmodel: "
198 # 200 #
199 # Run 1 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 caldb="irf" \ 205 caldb="irf" \
204 irf="cta_dummy_irf" \ 206 irf="cta_dummy_irf" \
205 - srcmdl="data/crab.xml" 207 + inmodel="data/crab.xml"
206 $ECHO -n "." 208 $ECHO -n "."
207 if [ -s "modmap1.fits" ] 209 if [ -s "modmap1.fits" ]
208 then 210 then
@@ -213,14 +215,15 @@ else @@ -213,14 +215,15 @@ else
213 fi 215 fi
214 # 216 #
215 # Run 3 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 caldb="irf" \ 221 caldb="irf" \
220 irf="cta_dummy_irf" \ 222 irf="cta_dummy_irf" \
221 - srcmdl="data/crab.xml" \ 223 + inmodel="data/crab.xml" \
222 ra=83.63 \ 224 ra=83.63 \
223 dec=22.01 \ 225 dec=22.01 \
  226 + rad=10.0 \
224 tmin=0.0 \ 227 tmin=0.0 \
225 tmax=1800.0 \ 228 tmax=1800.0 \
226 emin=0.1 \ 229 emin=0.1 \
@@ -243,12 +246,12 @@ else @@ -243,12 +246,12 @@ else
243 fi 246 fi
244 # 247 #
245 # Run 2 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 caldb="irf" \ 252 caldb="irf" \
250 irf="cta_dummy_irf" \ 253 irf="cta_dummy_irf" \
251 - srcmdl="data/model_background_cube.xml" 254 + inmodel="data/model_background_cube.xml"
252 $ECHO -n "." 255 $ECHO -n "."
253 if [ -s "modmap2.fits" ] 256 if [ -s "modmap2.fits" ]
254 then 257 then
@@ -264,8 +267,8 @@ $ECHO &quot; ok&quot; @@ -264,8 +267,8 @@ $ECHO &quot; ok&quot;
264 # Test ctselect 267 # Test ctselect
265 # ============= 268 # =============
266 $ECHO -n "Test ctselect: " 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 ra=83.63 \ 272 ra=83.63 \
270 dec=22.01 \ 273 dec=22.01 \
271 rad=10.0 \ 274 rad=10.0 \
@@ -290,9 +293,9 @@ $ECHO &quot; ok&quot; @@ -290,9 +293,9 @@ $ECHO &quot; ok&quot;
290 $ECHO -n "Test ctlike: " 293 $ECHO -n "Test ctlike: "
291 # 294 #
292 # Run 1 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 caldb="irf" \ 299 caldb="irf" \
297 irf="cta_dummy_irf" 300 irf="cta_dummy_irf"
298 $ECHO -n "." 301 $ECHO -n "."
@@ -305,9 +308,9 @@ else @@ -305,9 +308,9 @@ else
305 fi 308 fi
306 # 309 #
307 # Run 2 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 caldb="irf" \ 314 caldb="irf" \
312 irf="cta_dummy_irf" 315 irf="cta_dummy_irf"
313 $ECHO -n "." 316 $ECHO -n "."
@@ -320,9 +323,9 @@ else @@ -320,9 +323,9 @@ else
320 fi 323 fi
321 # 324 #
322 # Run 3 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 caldb="irf" \ 329 caldb="irf" \
327 irf="cta_dummy_irf" 330 irf="cta_dummy_irf"
328 $ECHO -n "." 331 $ECHO -n "."
@@ -335,9 +338,9 @@ else @@ -335,9 +338,9 @@ else
335 fi 338 fi
336 # 339 #
337 # Run 4 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 caldb="irf" \ 344 caldb="irf" \
342 irf="cta_dummy_irf" 345 irf="cta_dummy_irf"
343 $ECHO -n "." 346 $ECHO -n "."
@@ -355,12 +358,12 @@ $ECHO &quot; ok&quot; @@ -355,12 +358,12 @@ $ECHO &quot; ok&quot;
355 # Test cttsmap 358 # Test cttsmap
356 # ============ 359 # ============
357 $ECHO -n "Test cttsmap: " 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 srcname="Crab" \ 363 srcname="Crab" \
361 caldb="irf" \ 364 caldb="irf" \
362 irf="cta_dummy_irf" \ 365 irf="cta_dummy_irf" \
363 - outfile="tsmap.fits" \ 366 + outmap="tsmap.fits" \
364 nxpix=5 \ 367 nxpix=5 \
365 nypix=5 \ 368 nypix=5 \
366 binsz=0.05 \ 369 binsz=0.05 \
@@ -385,9 +388,9 @@ $ECHO &quot; ok&quot; @@ -385,9 +388,9 @@ $ECHO &quot; ok&quot;
385 $ECHO -n "Test ctexpcube: " 388 $ECHO -n "Test ctexpcube: "
386 # 389 #
387 # Run 1 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 caldb="irf" \ 394 caldb="irf" \
392 irf="cta_dummy_irf" \ 395 irf="cta_dummy_irf" \
393 emin=0.1 \ 396 emin=0.1 \
@@ -410,9 +413,9 @@ else @@ -410,9 +413,9 @@ else
410 fi 413 fi
411 # 414 #
412 # Run 2 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 caldb="irf" \ 419 caldb="irf" \
417 irf="cta_dummy_irf" 420 irf="cta_dummy_irf"
418 $ECHO -n "." 421 $ECHO -n "."
@@ -432,9 +435,9 @@ $ECHO &quot; ok&quot; @@ -432,9 +435,9 @@ $ECHO &quot; ok&quot;
432 $ECHO -n "Test ctpsfcube: " 435 $ECHO -n "Test ctpsfcube: "
433 # 436 #
434 # Run 1 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 caldb="irf" \ 441 caldb="irf" \
439 irf="cta_dummy_irf" \ 442 irf="cta_dummy_irf" \
440 emin=0.1 \ 443 emin=0.1 \
@@ -459,9 +462,9 @@ else @@ -459,9 +462,9 @@ else
459 fi 462 fi
460 # 463 #
461 # Run 2 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 caldb="irf" \ 468 caldb="irf" \
466 irf="cta_dummy_irf" \ 469 irf="cta_dummy_irf" \
467 amax=0.3 \ 470 amax=0.3 \
@@ -483,11 +486,13 @@ $ECHO &quot; ok&quot; @@ -483,11 +486,13 @@ $ECHO &quot; ok&quot;
483 $ECHO -n "Test ctbkgcube: " 486 $ECHO -n "Test ctbkgcube: "
484 # 487 #
485 # Run 1 488 # Run 1
486 -$ctbkgcube infile="data/crab_events.fits" \ 489 +$ctbkgcube inobs="data/crab_events.fits" \
487 inmodel="data/crab.xml" \ 490 inmodel="data/crab.xml" \
488 - cntmap="NONE" \  
489 - outfile="bkgcube1.fits" \ 491 + incube="NONE" \
  492 + outcube="bkgcube1.fits" \
490 outmodel="bkgcube1.xml" \ 493 outmodel="bkgcube1.xml" \
  494 + caldb="irf" \
  495 + irf="cta_dummy_irf" \
491 emin=0.1 \ 496 emin=0.1 \
492 emax=100.0 \ 497 emax=100.0 \
493 enumbins=20 \ 498 enumbins=20 \
@@ -515,10 +520,12 @@ else @@ -515,10 +520,12 @@ else
515 fi 520 fi
516 # 521 #
517 # Run 2 522 # Run 2
518 -$ctbkgcube infile="data/crab_events.fits" \ 523 +$ctbkgcube inobs="data/crab_events.fits" \
519 inmodel="data/crab.xml" \ 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 outmodel="bkgcube2.xml" 529 outmodel="bkgcube2.xml"
523 $ECHO -n "." 530 $ECHO -n "."
524 if [ -s "bkgcube2.fits" ] 531 if [ -s "bkgcube2.fits" ]
@@ -542,9 +549,9 @@ $ECHO &quot; ok&quot; @@ -542,9 +549,9 @@ $ECHO &quot; ok&quot;
542 # Test ctcubemask 549 # Test ctcubemask
543 # =============== 550 # ===============
544 $ECHO -n "Test ctcubemask: " 551 $ECHO -n "Test ctcubemask: "
545 -$ctcubemask infile="data/crab_cntmap.fits" \ 552 +$ctcubemask inobs="data/crab_cntmap.fits" \
546 regfile="data/exclusion.reg" \ 553 regfile="data/exclusion.reg" \
547 - outfile="filtered_cube.fits" \ 554 + outcube="filtered_cube.fits" \
548 ra=83.63 \ 555 ra=83.63 \
549 dec=22.01 \ 556 dec=22.01 \
550 rad=2.0 \ 557 rad=2.0 \
@@ -565,8 +572,8 @@ $ECHO &quot; ok&quot; @@ -565,8 +572,8 @@ $ECHO &quot; ok&quot;
565 # Test ctbutterfly 572 # Test ctbutterfly
566 # =============== 573 # ===============
567 $ECHO -n "Test ctbutterfly: " 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 srcname="Crab" \ 577 srcname="Crab" \
571 outfile="butterfly.txt" \ 578 outfile="butterfly.txt" \
572 caldb="irf" \ 579 caldb="irf" \
test/test_ctpsfcube.py
@@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite): @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite):
66 """ 66 """
67 # Set-up ctpsfcube 67 # Set-up ctpsfcube
68 psfcube = ctools.ctpsfcube() 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 psfcube["caldb"].string(self.caldb) 72 psfcube["caldb"].string(self.caldb)
73 psfcube["irf"].string(self.irf) 73 psfcube["irf"].string(self.irf)
74 psfcube["ebinalg"].string("LOG") 74 psfcube["ebinalg"].string("LOG")
test/test_ctselect.py
@@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 """ 64 """
65 # Set-up ctselect 65 # Set-up ctselect
66 select = ctools.ctselect() 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 select["ra"].real(83.63) 69 select["ra"].real(83.63)
70 select["dec"].real(22.01) 70 select["dec"].real(22.01)
71 select["rad"].real(3.0) 71 select["rad"].real(3.0)
test/test_ctskymap.py
@@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite):
64 """ 64 """
65 # Set-up ctskymap 65 # Set-up ctskymap
66 skymap = ctools.ctskymap() 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 skymap["emin"].real(0.1) 69 skymap["emin"].real(0.1)
70 skymap["emax"].real(100.0) 70 skymap["emax"].real(100.0)
71 skymap["nxpix"].integer(200) 71 skymap["nxpix"].integer(200)
test/test_cttsmap.py
@@ -67,12 +67,12 @@ class Test(gammalib.GPythonTestSuite): @@ -67,12 +67,12 @@ class Test(gammalib.GPythonTestSuite):
67 """ 67 """
68 # Set-up cttsmap 68 # Set-up cttsmap
69 tsmap = ctools.cttsmap() 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 tsmap["srcname"].string("Crab") 72 tsmap["srcname"].string("Crab")
73 tsmap["caldb"].string(self.caldb) 73 tsmap["caldb"].string(self.caldb)
74 tsmap["irf"].string(self.irf) 74 tsmap["irf"].string(self.irf)
75 - tsmap["outfile"].filename("tsmap.fits") 75 + tsmap["outmap"].filename("tsmap.fits")
76 tsmap["nxpix"].integer(5) 76 tsmap["nxpix"].integer(5)
77 tsmap["nypix"].integer(5) 77 tsmap["nypix"].integer(5)
78 tsmap["binsz"].real(0.02) 78 tsmap["binsz"].real(0.02)
test/test_pipelines.py
@@ -78,8 +78,8 @@ class Test(gammalib.GPythonTestSuite): @@ -78,8 +78,8 @@ class Test(gammalib.GPythonTestSuite):
78 78
79 # Simulate events 79 # Simulate events
80 sim = ctools.ctobssim() 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 sim["caldb"].string(caldb) 83 sim["caldb"].string(caldb)
84 sim["irf"].string(irf) 84 sim["irf"].string(irf)
85 sim["ra"].real(ra) 85 sim["ra"].real(ra)
@@ -98,8 +98,8 @@ class Test(gammalib.GPythonTestSuite): @@ -98,8 +98,8 @@ class Test(gammalib.GPythonTestSuite):
98 98
99 # Select events 99 # Select events
100 select = ctools.ctselect() 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 select["ra"].real(ra) 103 select["ra"].real(ra)
104 select["dec"].real(dec) 104 select["dec"].real(dec)
105 select["rad"].real(rad_select) 105 select["rad"].real(rad_select)
@@ -116,9 +116,9 @@ class Test(gammalib.GPythonTestSuite): @@ -116,9 +116,9 @@ class Test(gammalib.GPythonTestSuite):
116 116
117 # Perform maximum likelihood fitting 117 # Perform maximum likelihood fitting
118 like = ctools.ctlike() 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 like["caldb"].string(caldb) 122 like["caldb"].string(caldb)
123 like["irf"].string(irf) 123 like["irf"].string(irf)
124 self.test_try("Execute ctlike") 124 self.test_try("Execute ctlike")
@@ -149,7 +149,7 @@ class Test(gammalib.GPythonTestSuite): @@ -149,7 +149,7 @@ class Test(gammalib.GPythonTestSuite):
149 149
150 # Simulate events 150 # Simulate events
151 sim = ctools.ctobssim() 151 sim = ctools.ctobssim()
152 - sim["infile"].filename(model_name) 152 + sim["inmodel"].filename(model_name)
153 sim["caldb"].string(caldb) 153 sim["caldb"].string(caldb)
154 sim["irf"].string(irf) 154 sim["irf"].string(irf)
155 sim["ra"].real(ra) 155 sim["ra"].real(ra)
test/test_python.py
@@ -102,6 +102,7 @@ if __name__ == &#39;__main__&#39;: @@ -102,6 +102,7 @@ if __name__ == &#39;__main__&#39;:
102 102
103 # Run test suite 103 # Run test suite
104 success = suites.run() 104 success = suites.run()
  105 +
105 106
106 # Save test results 107 # Save test results
107 suites.save("reports/ctools.xml") 108 suites.save("reports/ctools.xml")