Mathieu Servillat / ctools

Commit 31f23bc7b8017fa83c4232cf528a79b7d671ee16

Authored by Jürgen Knödlseder
2 parents 430c2b2c 9119dabf

Merge branch 'integration' into devel

Inline Side-by-side
Showing 43 changed files with 1872 additions and 1682 deletions
cscripts/Makefile.am
  View file @31f23bc
... ... @@ -11,9 +11,10 @@ pythonscripts = __init__.py \
11 11  
12 12 # Parfiles to be distributed in $(prefix)/syspfiles
13 13 parfiles = cspull.par cssens.par cstsdist.par cscaldb.par csresmap.par \
14   - csobsdef.par csspec.par csobsinfo.par \
  14 + csobsdef.par csspec.par cslightcrv.par csobsinfo.par \
15 15 csmodelinfo.par csmodelmerge.par cstsmapmerge.par \
16   - csobs2caldb.par
  16 + csobs2caldb.par \
  17 + csfindobs.par csiactcopy.par csiactdata.par csiactobs.par
17 18  
18 19 # Scripts for which symbolic links will be created in bindir
19 20 linkscripts = cspull cssens cstsdist cscaldb csresmap \
... ...
cscripts/cscaldb.py
  View file @31f23bc
... ... @@ -37,7 +37,7 @@ class cscaldb(ctools.cscript):
37 37 into the "cscaldb.log" log file.
38 38 """
39 39  
40   - # Constructors and destructors
  40 + # Constructor
41 41 def __init__(self, *argv):
42 42 """
43 43 Constructor.
... ... @@ -46,28 +46,13 @@ class cscaldb(ctools.cscript):
46 46 self._name = "cscaldb"
47 47 self._version = "1.1.0"
48 48  
49   - # Initialise application
50   - if len(argv) == 0:
51   - ctools.cscript.__init__(self, self._name, self._version)
52   - elif len(argv) ==1:
53   - ctools.cscript.__init__(self, self._name, self._version, *argv)
54   - else:
55   - raise TypeError("Invalid number of arguments given.")
56   -
57   - # Set logger properties
58   - self._log_header()
59   - self._log.date(True)
  49 + # Initialise application by calling the appropriate class
  50 + # constructor.
  51 + self._init_cscript(argv)
60 52  
61 53 # Return
62 54 return
63 55  
64   - def __del__(self):
65   - """
66   - Destructor.
67   - """
68   - # Return
69   - return
70   -
71 56  
72 57 # Private methods
73 58 def _get_parameters(self):
... ... @@ -164,16 +149,6 @@ class cscaldb(ctools.cscript):
164 149  
165 150  
166 151 # Public methods
167   - def execute(self):
168   - """
169   - Execute the script.
170   - """
171   - # Run the script
172   - self.run()
173   -
174   - # Return
175   - return
176   -
177 152 def run(self):
178 153 """
179 154 Run the script.
... ... @@ -238,6 +213,19 @@ class cscaldb(ctools.cscript):
238 213 # Return
239 214 return
240 215  
  216 + def execute(self):
  217 + """
  218 + Execute the script.
  219 + """
  220 + # Open logfile
  221 + self._logFileOpen()
  222 +
  223 + # Run the script
  224 + self.run()
  225 +
  226 + # Return
  227 + return
  228 +
241 229  
242 230 # ======================== #
243 231 # Main routine entry point #
... ... @@ -247,8 +235,5 @@ if __name__ == '__main__':
247 235 # Create instance of application
248 236 app = cscaldb(sys.argv)
249 237  
250   - # Open logfile
251   - app._logFileOpen()
252   -
253 238 # Execute application
254 239 app.execute()
... ...
cscripts/csfindobs.par 0 → 100644
  View file @31f23bc
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csfindobs.par #
  8 +# Version: 1.1.0 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the csfindobs script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input and output parameters
  19 +#============================
  20 +datapath, s, a, $VHEFITS,,, "Path were data are located"
  21 +prodname, s, a, prod-name,,, "Name of FITS production (Run csiactdata to view your options)"
  22 +outfile, f, a, runlist.lis,,, "Output runlist file"
  23 +
  24 +#
  25 +# Script parameters
  26 +#==================
  27 +ra, r, a, 83.63,,, "Right ascension of selection region centre (deg)"
  28 +dec, r, a, 22.01,,, "Declination of selection region centre (deg)"
  29 +rad, r, a, 2.5,,, "Search radius of selection region (deg)"
  30 +min_qual, i, h, 0,0|1|2,, "Minimum data quality (0=perfect, 1=ok, 2=bad)"
  31 +expression, s, h, NONE,,, "Additional expression"
  32 +master_indx, s, h, master.json,,, "Name of master index file"
  33 +
  34 +#
  35 +# Standard parameters
  36 +#====================
  37 +chatter, i, h, 2,0,4, "Chattiness of output"
  38 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  39 +debug, b, h, no,,, "Debugging mode activated"
  40 +mode, s, h, ql,,, "Mode of automatic parameters"
  41 +logfile, f, h, csfindobs.log,,, "Log filename"
  42 +
  43 +
  44 +
... ...
cscripts/csfindobs.py
  View file @31f23bc
... ... @@ -33,7 +33,7 @@ class csfindobs(ctools.cscript):
33 33 Find observations from an IACT data store.
34 34 """
35 35  
36   - # Constructors and destructors
  36 + # Constructor
37 37 def __init__(self, *argv):
38 38 """
39 39 Constructor.
... ... @@ -44,68 +44,15 @@ class csfindobs(ctools.cscript):
44 44 self._verbose = False
45 45 self._datapath = os.getenv("VHEFITS","")
46 46  
47   - # Make sure that parfile exists
48   - self._parfile()
49   -
50   - # Initialise application
51   - if len(argv) == 0:
52   - ctools.cscript.__init__(self, self._name, self._version)
53   - elif len(argv) ==1:
54   - ctools.cscript.__init__(self, self._name, self._version, *argv)
55   - else:
56   - raise TypeError("Invalid number of arguments given.")
57   -
58   - # Set logger properties
59   - self._log_header()
60   - self._log.date(True)
  47 + # Initialise application by calling the appropriate class
  48 + # constructor.
  49 + self._init_cscript(argv)
61 50  
62 51 # Return
63 52 return
64 53  
65   - def __del__(self):
66   - """
67   - Destructor.
68   - """
69   - # Write separator into logger
70   - if self._logTerse():
71   - self._log("\n")
72   -
73   - # Return
74   - return
75 54  
76   - def _parfile(self):
77   - """
78   - Check if parfile exists. If parfile does not exist then create a
79   - default parfile. This kluge avoids shipping the cscript with a
80   - parfile.
81   - """
82   - # Set parfile name
83   - parfile = self._name+".par"
84   -
85   - try:
86   - pars = gammalib.GApplicationPars(parfile)
87   - except:
88   - # Signal if parfile was not found
89   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
90   -
91   - # Create default parfile
92   - pars = gammalib.GApplicationPars()
93   - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data are located"))
94   - pars.append(gammalib.GApplicationPar("prodname","s","a","prod-name","","","Name of FITS production (Run csiactdata to view your options)"))
95   - pars.append(gammalib.GApplicationPar("outfile","f","a","runlist.lis","","","Runlist outfile"))
96   - pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","","","Right ascension"))
97   - pars.append(gammalib.GApplicationPar("dec","r","a","22.01","","","Declination"))
98   - pars.append(gammalib.GApplicationPar("rad","r","a","2.5","","","Search radius"))
99   - pars.append(gammalib.GApplicationPar("min_qual","i","h","0","0|1|2","","Minimum data quality (0=perfect, 1=ok, 2=bad)"))
100   - pars.append(gammalib.GApplicationPar("expression","s","h","NONE","","","Additional expression"))
101   - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file"))
102   - pars.append_standard()
103   - pars.append(gammalib.GApplicationPar("logfile","f","h","csfindobs.log","","","Log filename"))
104   - pars.save(parfile)
105   -
106   - # Return
107   - return
108   -
  55 + # Private methods
109 56 def _get_parameters(self):
110 57 """
111 58 Get parameters from parfile and setup the observation.
... ... @@ -209,27 +156,7 @@ class csfindobs(ctools.cscript):
209 156 # Return
210 157 return
211 158  
212   - def obs_ids(self):
213   - """
214   - Return OBS IDs
215   - """
216   -
217   - # Return
218   - return self.runs
219   -
220   - def execute(self):
221   - """
222   - Execute the script.
223   - """
224   - # Run the script
225   - self.run()
226   -
227   - # Save residual map
228   - self.save(self._outfile)
229   -
230   - # Return
231   - return
232   -
  159 + # Public methods
233 160 def run(self):
234 161 """
235 162 Run the script.
... ... @@ -304,6 +231,22 @@ class csfindobs(ctools.cscript):
304 231 # Return
305 232 return
306 233  
  234 + def execute(self):
  235 + """
  236 + Execute the script.
  237 + """
  238 + # Open logfile
  239 + self._logFileOpen()
  240 +
  241 + # Run the script
  242 + self.run()
  243 +
  244 + # Save residual map
  245 + self.save(self._outfile)
  246 +
  247 + # Return
  248 + return
  249 +
307 250 def save(self, outfile):
308 251 """
309 252 Save.
... ... @@ -323,6 +266,14 @@ class csfindobs(ctools.cscript):
323 266 # Return
324 267 return
325 268  
  269 + def obs_ids(self):
  270 + """
  271 + Return OBS IDs
  272 + """
  273 +
  274 + # Return
  275 + return self._runs
  276 +
326 277  
327 278 # ======================== #
328 279 # Main routine entry point #
... ... @@ -332,8 +283,5 @@ if __name__ == '__main__':
332 283 # Create instance of application
333 284 app = csfindobs(sys.argv)
334 285  
335   - # Open logfile
336   - app._logFileOpen()
337   -
338 286 # Execute application
339 287 app.execute()
... ...
cscripts/csiactcopy.par 0 → 100644
  View file @31f23bc
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csiactcopy.par #
  8 +# Version: 1.1.0 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the csiactcopy script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input and output parameters
  19 +#============================
  20 +remote_master, f, a, /Volumes/mountpoint/master.json,,, "Location of remote master file"
  21 +prodname, s, a, iact-fits,,, "Name of FITS production to download"
  22 +runlist, f, h, NONE,,, "Input list of observation IDs"
  23 +outpath, s, a, /path/to/fits,,, "Destination path of FITS data"
  24 +
  25 +#
  26 +# Standard parameters
  27 +#====================
  28 +chatter, i, h, 2,0,4, "Chattiness of output"
  29 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  30 +debug, b, h, no,,, "Debugging mode activated"
  31 +mode, s, h, ql,,, "Mode of automatic parameters"
  32 +logfile, f, h, csiactcopy.log,,, "Log filename"
... ...
cscripts/csiactcopy.py
  View file @31f23bc
... ... @@ -34,6 +34,8 @@ class csiactcopy(ctools.cscript):
34 34 It can take a list of observation IDs to allow the download specific
35 35 observations. Index files get merged and updated accordingly.
36 36 """
  37 +
  38 + # Constructor
37 39 def __init__(self, *argv):
38 40 """
39 41 Constructor.
... ... @@ -43,60 +45,15 @@ class csiactcopy(ctools.cscript):
43 45 self._version = "1.1.0"
44 46 self._datapath = os.getenv("VHEFITS","")
45 47  
46   - # Make sure that parfile exists
47   - self._parfile()
48   -
49   - # Initialise application
50   - if len(argv) == 0:
51   - ctools.cscript.__init__(self, self._name, self._version)
52   - elif len(argv) == 1:
53   - ctools.cscript.__init__(self, self._name, self._version, *argv)
54   - else:
55   - raise TypeError("Invalid number of arguments given.")
56   -
57   - # Set logger properties
58   - self._log_header()
59   - self._log.date(True)
  48 + # Initialise application by calling the appropriate class
  49 + # constructor.
  50 + self._init_cscript(argv)
60 51  
61 52 # Return
62 53 return
63 54  
64   - def __del__(self):
65   - """
66   - Destructor.
67   - """
68   -
69   - # Return
70   - return
71   -
72   - def _parfile(self):
73   - """
74   - Check if parfile exists. If parfile does not exist then create a
75   - default parfile. This kluge avoids shipping the cscript with a
76   - parfile.
77   - """
78   - # Set parfile name
79   - parfile = self._name+".par"
80   -
81   - try:
82   - pars = gammalib.GApplicationPars(parfile)
83   - except:
84   - # Signal if parfile was not found
85   - print("Parfile "+parfile+" not found. Create default parfile.")
86   -
87   - # Create default parfile
88   - pars = gammalib.GApplicationPars()
89   - pars.append(gammalib.GApplicationPar("remote_master","f","a","/Volumes/mountpoint/master.json","","","Location of remote master file"))
90   - pars.append(gammalib.GApplicationPar("prodname","s","a","iact-fits","","","Name of FITS production to download"))
91   - pars.append(gammalib.GApplicationPar("outpath","s","a","/path/to/fits","","","Destination path of FITS data"))
92   - pars.append(gammalib.GApplicationPar("runlist","f","h","NONE","","","List of observation IDs"))
93   - pars.append_standard()
94   - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactcopy.log","","","Log filename"))
95   - pars.save(parfile)
96   -
97   - # Return
98   - return
99 55  
  56 + # Private methods
100 57 def _get_parameters(self):
101 58 """
102 59 Get parameters from parfile and setup the observation.
... ... @@ -114,20 +71,233 @@ class csiactcopy(ctools.cscript):
114 71 self._prodname = self["prodname"].string()
115 72 self._outpath = gammalib.expand_env(self["outpath"].string())
116 73 self._runlist = self["runlist"].filename()
117   -
  74 +
  75 + # Write input parameters into logger
  76 + if self._logTerse():
  77 + self._log_parameters()
  78 + self._log("\n")
  79 +
118 80 # Return
119 81 return
120 82  
121   - def execute(self):
  83 + def copy(self, source, clobber):
122 84 """
123   - Execute the script.
  85 + copy file to _outpath
124 86 """
125   - # Run the script
126   - self.run()
  87 + # Get file destination
  88 + destination = os.path.join(self._outpath,
  89 + os.path.relpath(source, self._remote_base))
  90 +
  91 + # Initialise return value
  92 + filesize = 0.0
  93 +
  94 + # Logging
  95 + if self._logVerbose():
  96 + self._log("\n")
  97 + self._log.header3("Copy file")
  98 + self._log.parformat("Source")
  99 + self._log(str(source))
  100 + self._log("\n")
  101 + self._log.parformat("Destination")
  102 + self._log(str(destination))
  103 + self._log("\n")
  104 + self._log.parformat("Already exists")
  105 + self._log(str(os.path.isfile(destination)))
  106 + self._log("\n")
  107 + self._log.parformat("Overwrite")
  108 + self._log(str(clobber))
  109 + self._log("\n")
  110 + self._log.parformat("Copying")
  111 +
  112 + # Flag if file destination is already available
  113 + is_file = os.path.isfile(destination)
  114 +
  115 + # check if file could be skipped because clobber=no
  116 + if is_file and clobber == False:
  117 + if self._logVerbose():
  118 + self._log("Skip (clobber=no)")
  119 + self._log("\n")
  120 +
  121 + # check if file could be skipped because it is the same file
  122 + elif is_file and os.path.samefile(destination, source):
  123 + if self._logVerbose():
  124 + self._log("Skip (same file)")
  125 + self._log("\n")
  126 +
  127 + else:
  128 + # Create directory if not existent
  129 + dest_dir = os.path.dirname(destination)
  130 + if not os.path.isdir(dest_dir):
  131 + os.makedirs(dest_dir)
  132 +
  133 + # Copy file (if source not identical to destination)
  134 + shutil.copy2(source, destination)
  135 +
  136 + # Get filesize
  137 + filesize = os.stat(destination).st_size
  138 +
  139 + # Logging
  140 + if self._logVerbose():
  141 + self._log("Done!")
  142 + self._log("\n")
  143 +
  144 + # Return
  145 + return filesize
  146 +
  147 +
  148 + def merge(self, localfits, remotefits, hduname, clobber):
  149 + """
  150 + merge remote and local fits file
  151 + If local fits file not present, a new one is created
  152 + """
  153 +
  154 + if self._logTerse():
  155 + self._log.parformat("Local file")
  156 + self._log(str(localfits))
  157 + self._log("\n")
  158 + self._log.parformat("Remote file")
  159 + self._log(str(remotefits))
  160 + self._log("\n")
  161 + self._log.parformat("HDU name")
  162 + self._log(str(hduname))
  163 + self._log("\n")
  164 + self._log.parformat("Conflict preference")
  165 + if clobber:
  166 + self._log("Remote")
  167 + else:
  168 + self._log("Local")
  169 + self._log("\n")
  170 +
  171 + # Check if local fits file is available
  172 + if os.path.isfile(localfits):
  173 + local = gammalib.GFits(str(localfits))
  174 + local_hdu = local[hduname]
  175 + else:
  176 + # Otherwise load remote fits file and delete rows
  177 + local = gammalib.GFits(str(remotefits))
  178 + local_hdu = local[hduname]
  179 + local_hdu.remove_rows(0,local_hdu.nrows())
  180 +
  181 + # find local obs_id
  182 + local_obs = []
  183 + lobs_col = local_hdu["OBS_ID"]
  184 + for obs_id in lobs_col:
  185 + local_obs.append(obs_id)
  186 +
  187 + # Create selection expression for opening remote fits file
  188 + selection = ""
  189 + if len(self._runs):
  190 + selection += "["+hduname+"]["
  191 + for run in self._runs:
  192 + selection += "OBS_ID=="+str(run)
  193 + selection += "||"
  194 + selection = selection[:-2]
  195 + selection += "]"
  196 +
  197 + # Create file name
  198 + remotefile = remotefits + selection
  199 +
  200 + # open remote fits file
  201 + remote = gammalib.GFits(str(remotefile))
  202 + remote_hdu = remote[hduname]
  203 +
  204 + # Get remote obs_id
  205 + remote_obs = []
  206 + robs_col = remote_hdu["OBS_ID"]
  207 + for obs_id in robs_col:
  208 + remote_obs.append(obs_id)
  209 +
  210 + if self._logTerse():
  211 + self._log.parformat("Remote entries")
  212 + self._log(str(len(remote_obs)))
  213 + self._log("\n")
  214 + self._log.parformat("Local entries")
  215 + self._log(str(len(local_obs)))
  216 + self._log("\n")
  217 +
  218 + # initialise counter for logging
  219 + removed_rows = 0
  220 +
  221 + # If clobber=yes, we use the remote content to overwrite
  222 + # corresponding local content
  223 + if clobber:
  224 +
  225 + # Loop over remote obs_ids
  226 + for remote_obs_id in remote_obs:
  227 +
  228 + # Check if remote obs_id is present locally
  229 + if remote_obs_id in local_obs:
  230 +
  231 + # Remove local obs_id entries
  232 + table_has_obsid = True
  233 + while table_has_obsid:
  234 + for i in range(local_hdu.nrows()):
  235 + if remote_obs_id == local_hdu["OBS_ID"][i]:
  236 + local_hdu.remove_rows(i,1)
  237 + removed_rows += 1
  238 + break
  239 + table_has_obsid = False
  240 + for i in range(local_hdu.nrows()):
  241 + if remote_obs_id == local_hdu["OBS_ID"][i]:
  242 + table_has_obsid = True
  243 + break
  244 +
  245 + # If clobber=false, we keep the local content unchanged and just
  246 + # expand by remote content
  247 + else:
  248 + # Loop over remote obs_ids
  249 + for local_obs_id in local_obs:
  250 +
  251 + # Check if local obs_id is present remotely
  252 + if local_obs_id in remote_obs:
  253 +
  254 + # Remove remote obs_id entries
  255 + table_has_obsid = True
  256 + while table_has_obsid:
  257 + for i in range(remote_hdu.nrows()):
  258 + if local_obs_id == remote_hdu["OBS_ID"][i]:
  259 + remote_hdu.remove_rows(i,1)
  260 + removed_rows += 1
  261 + break
  262 + table_has_obsid = False
  263 + for i in range(remote_hdu.nrows()):
  264 + if local_obs_id == remote_hdu["OBS_ID"][i]:
  265 + table_has_obsid = True
  266 + break
  267 +
  268 + # Store number of local rows
  269 + old_local_rows = local_hdu.nrows()
  270 +
  271 + # Add rmeote HDUs
  272 + local_hdu.insert_rows(old_local_rows, remote_hdu.nrows())
  273 +
  274 + if self._logTerse():
  275 + if clobber:
  276 + self._log.parformat("Removed local rows")
  277 + else:
  278 + self._log.parformat("Skipped remote rows")
  279 + self._log(str(removed_rows))
  280 + self._log("\n")
127 281  
  282 + # Loop over columns
  283 + for i in range(local_hdu.ncols()):
  284 +
  285 + # Get local and remote columns
  286 + local_col = local_hdu[i]
  287 + remote_col = remote_hdu[i]
  288 +
  289 + # Loop over entries and merge
  290 + for j in range(remote_col.length()):
  291 + local_col[j+old_local_rows] = remote_col[j]
  292 +
  293 + # Save local fits file
  294 + local.saveto(str(localfits), True)
  295 +
128 296 # Return
129 297 return
130   -
  298 +
  299 +
  300 + # Public methods
131 301 def run(self):
132 302 """
133 303 Run the script.
... ... @@ -138,11 +308,6 @@ class csiactcopy(ctools.cscript):
138 308 # Get parameters
139 309 self._get_parameters()
140 310  
141   - # Write input parameters into logger
142   - if self._logTerse():
143   - self._log_parameters()
144   - self._log("\n")
145   -
146 311 # Make destination directory if not available
147 312 if not os.path.isdir(self._outpath):
148 313 os.makedirs(self._outpath)
... ... @@ -516,221 +681,19 @@ class csiactcopy(ctools.cscript):
516 681 # Return
517 682 return
518 683  
519   - def copy(self, source, clobber):
520   - """
521   - copy file to m_outpath
  684 + def execute(self):
522 685 """
523   - # Get file destination
524   - destination = os.path.join(self._outpath,
525   - os.path.relpath(source, self._remote_base))
526   -
527   - # Initialise return value
528   - filesize = 0.0
529   -
530   - # Logging
531   - if self._logVerbose():
532   - self._log("\n")
533   - self._log.header3("Copy file")
534   - self._log.parformat("Source")
535   - self._log(str(source))
536   - self._log("\n")
537   - self._log.parformat("Destination")
538   - self._log(str(destination))
539   - self._log("\n")
540   - self._log.parformat("Already exists")
541   - self._log(str(os.path.isfile(destination)))
542   - self._log("\n")
543   - self._log.parformat("Overwrite")
544   - self._log(str(clobber))
545   - self._log("\n")
546   - self._log.parformat("Copying")
547   -
548   - # Flag if file destination is already available
549   - is_file = os.path.isfile(destination)
550   -
551   - # check if file could be skipped because clobber=no
552   - if is_file and clobber == False:
553   - if self._logVerbose():
554   - self._log("Skip (clobber=no)")
555   - self._log("\n")
556   -
557   - # check if file could be skipped because it is the same file
558   - elif is_file and os.path.samefile(destination, source):
559   - if self._logVerbose():
560   - self._log("Skip (same file)")
561   - self._log("\n")
562   -
563   - else:
564   - # Create directory if not existent
565   - dest_dir = os.path.dirname(destination)
566   - if not os.path.isdir(dest_dir):
567   - os.makedirs(dest_dir)
568   -
569   - # Copy file (if source not identical to destination)
570   - shutil.copy2(source, destination)
571   -
572   - # Get filesize
573   - filesize = os.stat(destination).st_size
574   -
575   - # Logging
576   - if self._logVerbose():
577   - self._log("Done!")
578   - self._log("\n")
579   -
580   - # Return
581   - return filesize
582   -
583   -
584   - def merge(self, localfits, remotefits, hduname, clobber):
  686 + Execute the script.
585 687 """
586   - merge remote and local fits file
587   - If local fits file not present, a new one is created
588   - """
589   -
590   - if self._logTerse():
591   - self._log.parformat("Local file")
592   - self._log(str(localfits))
593   - self._log("\n")
594   - self._log.parformat("Remote file")
595   - self._log(str(remotefits))
596   - self._log("\n")
597   - self._log.parformat("HDU name")
598   - self._log(str(hduname))
599   - self._log("\n")
600   - self._log.parformat("Conflict preference")
601   - if clobber:
602   - self._log("Remote")
603   - else:
604   - self._log("Local")
605   - self._log("\n")
606   -
607   - # Check if local fits file is available
608   - if os.path.isfile(localfits):
609   - local = gammalib.GFits(str(localfits))
610   - local_hdu = local[hduname]
611   - else:
612   - # Otherwise load remote fits file and delete rows
613   - local = gammalib.GFits(str(remotefits))
614   - local_hdu = local[hduname]
615   - local_hdu.remove_rows(0,local_hdu.nrows())
616   -
617   - # find local obs_id
618   - local_obs = []
619   - lobs_col = local_hdu["OBS_ID"]
620   - for obs_id in lobs_col:
621   - local_obs.append(obs_id)
622   -
623   - # Create selection expression for opening remote fits file
624   - selection = ""
625   - if len(self._runs):
626   - selection += "["+hduname+"]["
627   - for run in self._runs:
628   - selection += "OBS_ID=="+str(run)
629   - selection += "||"
630   - selection = selection[:-2]
631   - selection += "]"
632   -
633   - # Create file name
634   - remotefile = remotefits + selection
635   -
636   - # open remote fits file
637   - remote = gammalib.GFits(str(remotefile))
638   - remote_hdu = remote[hduname]
639   -
640   - # Get remote obs_id
641   - remote_obs = []
642   - robs_col = remote_hdu["OBS_ID"]
643   - for obs_id in robs_col:
644   - remote_obs.append(obs_id)
645   -
646   - if self._logTerse():
647   - self._log.parformat("Remote entries")
648   - self._log(str(len(remote_obs)))
649   - self._log("\n")
650   - self._log.parformat("Local entries")
651   - self._log(str(len(local_obs)))
652   - self._log("\n")
653   -
654   - # initialise counter for logging
655   - removed_rows = 0
656   -
657   - # If clobber=yes, we use the remote content to overwrite
658   - # corresponding local content
659   - if clobber:
660   -
661   - # Loop over remote obs_ids
662   - for remote_obs_id in remote_obs:
663   -
664   - # Check if remote obs_id is present locally
665   - if remote_obs_id in local_obs:
666   -
667   - # Remove local obs_id entries
668   - table_has_obsid = True
669   - while table_has_obsid:
670   - for i in range(local_hdu.nrows()):
671   - if remote_obs_id == local_hdu["OBS_ID"][i]:
672   - local_hdu.remove_rows(i,1)
673   - removed_rows += 1
674   - break
675   - table_has_obsid = False
676   - for i in range(local_hdu.nrows()):
677   - if remote_obs_id == local_hdu["OBS_ID"][i]:
678   - table_has_obsid = True
679   - break
680   -
681   - # If clobber=false, we keep the local content unchanged and just
682   - # expand by remote content
683   - else:
684   - # Loop over remote obs_ids
685   - for local_obs_id in local_obs:
686   -
687   - # Check if local obs_id is present remotely
688   - if local_obs_id in remote_obs:
689   -
690   - # Remove remote obs_id entries
691   - table_has_obsid = True
692   - while table_has_obsid:
693   - for i in range(remote_hdu.nrows()):
694   - if local_obs_id == remote_hdu["OBS_ID"][i]:
695   - remote_hdu.remove_rows(i,1)
696   - removed_rows += 1
697   - break
698   - table_has_obsid = False
699   - for i in range(remote_hdu.nrows()):
700   - if local_obs_id == remote_hdu["OBS_ID"][i]:
701   - table_has_obsid = True
702   - break
703   -
704   - # Store number of local rows
705   - old_local_rows = local_hdu.nrows()
706   -
707   - # Add rmeote HDUs
708   - local_hdu.insert_rows(old_local_rows, remote_hdu.nrows())
709   -
710   - if self._logTerse():
711   - if clobber:
712   - self._log.parformat("Removed local rows")
713   - else:
714   - self._log.parformat("Skipped remote rows")
715   - self._log(str(removed_rows))
716   - self._log("\n")
  688 + # Open logfile
  689 + self._logFileOpen()
  690 +
  691 + # Run the script
  692 + self.run()
717 693  
718   - # Loop over columns
719   - for i in range(local_hdu.ncols()):
720   -
721   - # Get local and remote columns
722   - local_col = local_hdu[i]
723   - remote_col = remote_hdu[i]
724   -
725   - # Loop over entries and merge
726   - for j in range(remote_col.length()):
727   - local_col[j+old_local_rows] = remote_col[j]
728   -
729   - # Save local fits file
730   - local.saveto(str(localfits), True)
731   -
732 694 # Return
733 695 return
  696 +
734 697  
735 698 # ======================== #
736 699 # Main routine entry point #
... ... @@ -740,8 +703,5 @@ if __name__ == '__main__':
740 703 # Create instance of application
741 704 app = csiactcopy(sys.argv)
742 705  
743   - # Open logfile
744   - app._logFileOpen()
745   -
746 706 # Execute application
747 707 app.execute()
... ...
cscripts/csiactdata.par 0 → 100644
  View file @31f23bc
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csiactdata.par #
  8 +# Version: 1.1.0 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the csiactdata script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input and output parameters
  19 +#============================
  20 +datapath, s, a, $VHEFITS,,, "Path were data is located"
  21 +master_indx, s, h, master.json,,, "Name of master index file"
  22 +
  23 +#
  24 +# Standard parameters
  25 +#====================
  26 +chatter, i, h, 2,0,4, "Chattiness of output"
  27 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  28 +debug, b, h, no,,, "Debugging mode activated"
  29 +mode, s, h, ql,,, "Mode of automatic parameters"
  30 +logfile, f, h, csiactdata.log,,, "Log filename"
  31 +
... ...
cscripts/csiactdata.py
  View file @31f23bc
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # Inspection of IACT data storage
4 4 #
5   -# Copyright (C) 2015 Michael Mayer
  5 +# Copyright (C) 2015-2016 Michael Mayer
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -35,6 +35,8 @@ class csiactdata(ctools.cscript):
35 35 file or prints it on the screen. These information can be used as input
36 36 for 'csiactobs'.
37 37 """
  38 +
  39 + # Constructor
38 40 def __init__(self, *argv):
39 41 """
40 42 Constructor.
... ... @@ -45,57 +47,15 @@ class csiactdata(ctools.cscript):
45 47 self._datapath = os.getenv("VHEFITS","")
46 48 self._prodnames = []
47 49  
48   - # Make sure that parfile exists
49   - file = self._parfile()
50   -
51   - # Initialise application
52   - if len(argv) == 0:
53   - ctools.cscript.__init__(self, self._name, self._version)
54   - elif len(argv) == 1:
55   - ctools.cscript.__init__(self, self._name, self._version, *argv)
56   - else:
57   - raise TypeError("Invalid number of arguments given.")
58   -
59   - # Set logger properties
60   - self._log_header()
61   - self._log.date(True)
  50 + # Initialise application by calling the appropriate class
  51 + # constructor.
  52 + self._init_cscript(argv)
62 53  
63 54 # Return
64 55 return
65 56  
66   - def __del__(self):
67   - """
68   - Destructor.
69   - """
70   - # Return
71   - return
72   -
73   - def _parfile(self):
74   - """
75   - Check if parfile exists. If parfile does not exist then create a
76   - default parfile. This kluge avoids shipping the cscript with a
77   - parfile.
78   - """
79   - # Set parfile name
80   - parfile = self._name+".par"
81   -
82   - try:
83   - pars = gammalib.GApplicationPars(parfile)
84   - except:
85   - # Signal if parfile was not found
86   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
87   -
88   - # Create default parfile
89   - pars = gammalib.GApplicationPars()
90   - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data is located"))
91   - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file"))
92   - pars.append_standard()
93   - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactdata.log","","","Log filename"))
94   - pars.save(parfile)
95   -
96   - # Return
97   - return
98 57  
  58 + # Private methods
99 59 def get_parameters(self):
100 60 """
101 61 Get parameters from parfile and setup the observation.
... ... @@ -119,16 +79,7 @@ class csiactdata(ctools.cscript):
119 79 # Return
120 80 return
121 81  
122   - def execute(self):
123   - """
124   - Execute the script.
125   - """
126   - # Run the script
127   - self.run()
128   -
129   - # Return
130   - return
131   -
  82 + # Public methods
132 83 def run(self):
133 84 """
134 85 Run the script.
... ... @@ -222,6 +173,20 @@ class csiactdata(ctools.cscript):
222 173 # Return
223 174 return
224 175  
  176 + def execute(self):
  177 + """
  178 + Execute the script.
  179 + """
  180 + # Open logfile
  181 + self._logFileOpen()
  182 +
  183 + # Run the script
  184 + self.run()
  185 +
  186 + # Return
  187 + return
  188 +
  189 +
225 190 def names(self):
226 191 """
227 192 Return available FITS production names
... ... @@ -239,8 +204,5 @@ if __name__ == '__main__':
239 204 # Create instance of application
240 205 app = csiactdata(sys.argv)
241 206  
242   - # Open logfile
243   - app._logFileOpen()
244   -
245 207 # Execute application
246 208 app.execute()
... ...
cscripts/csiactobs.par 0 → 100644
  View file @31f23bc
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csiactobs.par #
  8 +# Version: 1.1.0 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the csiactobs script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input and output parameters
  19 +#============================
  20 +datapath, s, a, $VHEFITS,,, "Path were data is located"
  21 +prodname, s, a, ,,, "Data storage name"
  22 +infile, f, a, runlist.lis,,, "Input runlist file"
  23 +inmodel, f, h, NONE,,, "Input model definition XML file (optional)"
  24 +master_indx, s, h, master.json,,, "Input master index file"
  25 +outobs, f, a, obs.xml,,, "Output observation definition XML file"
  26 +outmodel, f, a, bgmodels.xml,,, "Output model definition XML file"
  27 +
  28 +#
  29 +# Script parameters
  30 +#==================
  31 +bkgpars, i, a, 1,,, "Number of free parameters per background model"
  32 +bkg_scale, b, h, yes,,, "Apply background scaling factor from observation index file if available?"
  33 +ev_hiera, s, h, events,,, "Hierarchy of event formats"
  34 +aeff_hiera, s, h, aeff_2d,,, "Hierarchy of effective area formats"
  35 +psf_hiera, s, h, psf_king|psf_3gauss,,, "Hierarchy of psf formats"
  36 +edisp_hiera, s, h, edisp_2d,,, "Hierarchy of energy dispersion formats"
  37 +bkg_hiera, s, h, bkg_3d,,, "Hierarchy of background formats"
  38 +bkg_mod_hiera, s, h, irf|aeff|gauss,,, "Hierarchy of background models"
  39 +bkg_gauss_norm, r, h, 1e-8,,, "Input normalisation for Gaussian background"
  40 +bkg_gauss_index, r, h, -2.0,,, "Input spectral index for Gaussian background"
  41 +bkg_gauss_sigma, r, h, 2.5,,, "Input sigma for Gaussian background"
  42 +bkg_aeff_norm, r, h, 1e-14,,, "Input normalisation for effective area background"
  43 +bkg_aeff_index, r, h, -2.0,,, "Input spectral index for effective area background"
  44 +bkg_range_factor, r, h, 100.0,,, "Factor to determine range of background normalisation"
  45 +
  46 +#
  47 +# Standard parameters
  48 +#====================
  49 +chatter, i, h, 2,0,4, "Chattiness of output"
  50 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  51 +debug, b, h, no,,, "Debugging mode activated"
  52 +mode, s, h, ql,,, "Mode of automatic parameters"
  53 +logfile, f, h, csiactobs.log,,, "Log filename"
  54 +
... ...
cscripts/csiactobs.py
  View file @31f23bc
... ... @@ -37,6 +37,8 @@ class csiactobs(ctools.cscript):
37 37 described here:
38 38 http://gamma-astro-data-formats.readthedocs.org/en/latest/
39 39 """
  40 +
  41 + # Constructor
40 42 def __init__(self, *argv):
41 43 """
42 44 Constructor.
... ... @@ -53,75 +55,15 @@ class csiactobs(ctools.cscript):
53 55 self._models = gammalib.GModels()
54 56 self._runlist = []
55 57  
56   - # Make sure that parfile exists
57   - file = self._parfile()
58   -
59   - # Initialise application
60   - if len(argv) == 0:
61   - ctools.cscript.__init__(self, self._name, self._version)
62   - elif len(argv) == 1:
63   - ctools.cscript.__init__(self, self._name, self._version, *argv)
64   - else:
65   - raise TypeError("Invalid number of arguments given.")
66   -
67   - # Set logger properties
68   - self._log_header()
69   - self._log.date(True)
70   -
71   - # Return
72   - return
  58 + # Initialise application by calling the appropriate class
  59 + # constructor.
  60 + self._init_cscript(argv)
73 61  
74   - def __del__(self):
75   - """
76   - Destructor.
77   - """
78 62 # Return
79 63 return
80 64  
81   - def _parfile(self):
82   - """
83   - Check if parfile exists. If parfile does not exist then create a
84   - default parfile. This kluge avoids shipping the cscript with a parfile.
85   - """
86   - # Set parfile name
87   - parfile = self._name+".par"
88   -
89   - try:
90   - pars = gammalib.GApplicationPars(parfile)
91   - except:
92   - # Signal if parfile was not found
93   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
94   -
95   - # Create default parfile
96   - pars = gammalib.GApplicationPars()
97   - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data is located"))
98   - pars.append(gammalib.GApplicationPar("prodname","s","a","","","","Data storage name (Run csiactdata to view your options)"))
99   - pars.append(gammalib.GApplicationPar("infile","f","a","runlist.lis","","","Runlist file"))
100   - pars.append(gammalib.GApplicationPar("inmodel","f","h","NONE","","","Input model XML file (optional)"))
101   - pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile"))
102   - pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file"))
103   - pars.append(gammalib.GApplicationPar("bkgpars","i","a","1","","","Number of free parameters per background model"))
104   - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file"))
105   - pars.append(gammalib.GApplicationPar("bkg_scale","b","h","yes","","","Specifies whether the background scaling factor from the observation index file should be applied if available. "))
106   - pars.append(gammalib.GApplicationPar("ev_hiera","s","h","events","","","Hierarchy of event formats"))
107   - pars.append(gammalib.GApplicationPar("aeff_hiera","s","h","aeff_2d","","","Hierarchy of effective area formats"))
108   - pars.append(gammalib.GApplicationPar("psf_hiera","s","h","psf_king|psf_3gauss","","","Hierarchy of psf formats"))
109   - pars.append(gammalib.GApplicationPar("edisp_hiera","s","h","edisp_2d","","","Hierarchy of energy dispersion formats"))
110   - pars.append(gammalib.GApplicationPar("bkg_hiera","s","h","bkg_3d","","","Hierarchy of background formats"))
111   - pars.append(gammalib.GApplicationPar("bkg_mod_hiera","s","h","irf|aeff|gauss","","","Hierarchy of background models"))
112   - pars.append(gammalib.GApplicationPar("bkg_gauss_norm","r","h","1e-8","","","Input normalisation for Gaussian background"))
113   - pars.append(gammalib.GApplicationPar("bkg_gauss_index","r","h","-2.0","","","Input spectral index for Gaussian background"))
114   - pars.append(gammalib.GApplicationPar("bkg_gauss_sigma","r","h","2.5","","","Input sigma for Gaussian background"))
115   - pars.append(gammalib.GApplicationPar("bkg_aeff_norm","r","h","1e-14","","","Input normalisation for effective area background"))
116   - pars.append(gammalib.GApplicationPar("bkg_aeff_index","r","h","-2.0","","","Input spectral index for effective area background"))
117   - pars.append(gammalib.GApplicationPar("bkg_range_factor","r","h","100.0","","","Factor to determine range of background normalisation"))
118   - pars.append_standard()
119   - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactobs.log","","","Log filename"))
120   - pars.save(parfile)
121   -
122   - # Return
123   - return
124 65  
  66 + # Private methods
125 67 def _get_parameters(self):
126 68 """
127 69 Get parameters from parfile and setup the observation.
... ... @@ -366,43 +308,8 @@ class csiactobs(ctools.cscript):
366 308 # Return model
367 309 return model
368 310  
369   - def ebounds(self):
370   - """
371   - Returns runlist energy range
372   - """
373   - return self._ebounds
374   -
375   - def obs(self):
376   - """
377   - Returns GObservations object
378   - """
379   - # Initialise observations
380   - obs = gammalib.GObservations()
381   -
382   - # Read observations if XML is filled
383   - if self._xml.size():
384   - obs.read(self._xml)
385   -
386   - # Assign models
387   - obs.models(self._models)
388   -
389   - # Return observations
390   - return obs
391   -
392   -
393   - def execute(self):
394   - """
395   - Execute the script.
396   - """
397   - # Run the script
398   - self.run()
399   -
400   - # Save residual map
401   - self.save(self._outobs, self._clobber)
402   -
403   - # Return
404   - return
405 311  
  312 + # Public methods
406 313 def run(self):
407 314 """
408 315 Run the script.
... ... @@ -665,6 +572,22 @@ class csiactobs(ctools.cscript):
665 572 # Return
666 573 return
667 574  
  575 + def execute(self):
  576 + """
  577 + Execute the script.
  578 + """
  579 + # Open logfile
  580 + self._logFileOpen()
  581 +
  582 + # Run the script
  583 + self.run()
  584 +
  585 + # Save residual map
  586 + self.save(self._outobs, self._clobber)
  587 +
  588 + # Return
  589 + return
  590 +
668 591 def save(self,outfile,clobber):
669 592  
670 593 # Save observation XML file
... ... @@ -676,6 +599,29 @@ class csiactobs(ctools.cscript):
676 599 # Return
677 600 return
678 601  
  602 + def ebounds(self):
  603 + """
  604 + Returns runlist energy range
  605 + """
  606 + return self._ebounds
  607 +
  608 + def obs(self):
  609 + """
  610 + Returns GObservations object
  611 + """
  612 + # Initialise observations
  613 + obs = gammalib.GObservations()
  614 +
  615 + # Read observations if XML is filled
  616 + if self._xml.size():
  617 + obs.read(self._xml)
  618 +
  619 + # Assign models
  620 + obs.models(self._models)
  621 +
  622 + # Return observations
  623 + return obs
  624 +
679 625  
680 626 # ======================== #
681 627 # Main routine entry point #
... ... @@ -685,8 +631,5 @@ if __name__ == '__main__':
685 631 # Create instance of application
686 632 app = csiactobs(sys.argv)
687 633  
688   - # Open logfile
689   - app._logFileOpen()
690   -
691 634 # Execute application
692 635 app.execute()
... ...
cscripts/cslightcrv.par 0 → 100644
  View file @31f23bc
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: cslightcrv.par #
  8 +# Version: 1.1.0 #
  9 +# #
  10 +# Author: Juergen Knoedlseder #
  11 +# Jurgen.Knodlseder@irap.omp.eu #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the cslightcrv script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input and output parameters
  19 +#============================
  20 +inobs, f, h, events.fits,,, "Input event list, counts cube, or observation definition XML file"
  21 +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file"
  22 +srcname, s, a, Crab,,, "Source name"
  23 +expcube, f, a, NONE,,, "Input exposure cube file (only needed for stacked analysis)"
  24 +psfcube, f, a, NONE,,, "Input PSF cube file (only needed for stacked analysis)"
  25 +edispcube, f, a, NONE,,, "Input energy dispersion cube file (only needed for stacked analysis)"
  26 +bkgcube, f, a, NONE,,, "Input background cube file (only needed for stacked analysis)"
  27 +caldb, s, a, prod2,,, "Calibration database"
  28 +irf, s, a, South_0.5h,,, "Instrument response function"
  29 +edisp, b, h, no,,, "Apply energy dispersion?"
  30 +outfile, f, a, lightcurve.fits,,, "Output light curve file"
  31 +
  32 +#
  33 +# Script parameters
  34 +#==================
  35 +tbinalg, s, a, GTI,FILE|LIN|GTI,, "Algorithm for defining time bins"
  36 +tmin, r, a, 51544.5,,,"Lightcurve start time (MJD])"
  37 +tmax, r, a, 51544.6,,,"Lightcurve stop time (MJD)"
  38 +tbins, i, a, 5,1,10000, "Number of time bins"
  39 +tbinfile, f, a, NONE,,, "File defining the time binning"
  40 +emin, r, a, 0.1,,, "Lower energy limit of events (TeV)"
  41 +emax, r, a, 100.0,,, "Upper energy limit of events (TeV)"
  42 +enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)"
  43 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System"
  44 +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
  45 +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
  46 +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
  47 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  48 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  49 +binsz, r, a, 0.02,,, "Pixel size (deg/pixel)"
  50 +calc_ts, b, h, yes,,, "Compute TS value for each bin?"
  51 +calc_ulim, b, h, yes,,, "Compute upper limit for each bin?"
  52 +fix_srcs, b, h, yes,,, "Fix other sky model parameters?"
  53 +fix_bkg, b, h, no,,, "Fix background model parameters?"
  54 +
  55 +#
  56 +# Standard parameters
  57 +#====================
  58 +chatter, i, h, 2,0,4, "Chattiness of output"
  59 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  60 +debug, b, h, no,,, "Debugging mode activated"
  61 +mode, s, h, ql,,, "Mode of automatic parameters"
  62 +logfile, f, h, cslightcrv.log,,, "Log filename"
  63 +
  64 +
  65 +
... ...
cscripts/cslightcrv.py
  View file @31f23bc
... ... @@ -30,159 +30,70 @@ class cslightcrv(ctools.cscript):
30 30 """
31 31 Generates a lightcurve.
32 32  
33   - This class implements the creation of a light curve. It derives from
34   - the ctools.cscript class which provides support for parameter files,
35   - command line arguments, and logging. In that way the Python script
36   - behaves just as a regular ctool.
  33 + The cslightcrv class generates a light curve for Imaging Air Cherenkov
  34 + Telescope event data by performing a maximum likelihood fit using
  35 + ctlike in a series of time bins. The time bins can be either
  36 + specified in an ASCII file, as an interval divided into equally
  37 + sized time bins, or can be taken from the Good Time Intervals of the
  38 + observation(s).
  39 +
  40 + The format of the ASCII file is one row per time bin, each specifying
  41 + the start of stop value of the bin, separated by a whitespace. The
  42 + times are given in Modified Julian Days (MJD).
  43 +
  44 + Examples:
  45 + >>> lcrv = cslightcrv()
  46 + >>> lcrv.run()
  47 + >>> ... (querying for parameters) ...
  48 + >>> fits = lcrv.lightcurve()
  49 + Generates a light curve and retrieves the results in
  50 + a FITS file.
  51 +
  52 + >>> lcrv = cslightcrv()
  53 + >>> lcrv.execute()
  54 + >>> ... (querying for parameters) ...
  55 + Generates a light curve and saves results in a FITS file.
  56 +
  57 + >>> lcrv = cslightcrv(obs)
  58 + >>> lcrv.execute()
  59 + >>> ... (querying for parameters) ...
  60 + Generates a light curve from the observations in an
  61 + observation container and saves results in a FITS file.
37 62 """
38 63  
39   - # Constructors and destructors
  64 + # Constructor
40 65 def __init__(self, *argv):
41 66 """
42 67 Constructor.
43 68 """
44   -
45 69 # Set name
46 70 self._name = "cslightcrv"
47 71 self._version = "1.1.0"
48 72  
49 73 # Initialise some members
50   - self._outfile = ""
51   -
52   - # Initialise observation
53   - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
54   - self._obs = argv[0]
55   - argv = argv[1:]
56   - else:
57   - self._obs = gammalib.GObservations()
58   -
59   - # Make sure that parfile exists
60   - self._parfile()
61   -
62   - # Initialise application
63   - if len(argv) == 0:
64   - ctools.cscript.__init__(self, self._name, self._version)
65   - elif len(argv) ==1:
66   - ctools.cscript.__init__(self, self._name, self._version, *argv)
67   - else:
68   - raise TypeError("Invalid number of arguments given.")
  74 + self._srcname = ""
  75 + self._tbins = gammalib.GGti()
  76 + self._stacked = False
69 77  
70   - # Set logger properties
71   - self._log_header()
72   - self._log.date(True)
  78 + # Initialise observation container from constructor arguments.
  79 + self._obs = self._set_input_obs(argv)
73 80  
74   - # Return
75   - return
  81 + # Initialise application by calling the appropriate class
  82 + # constructor.
  83 + self._init_cscript(argv)
76 84  
77   - def __del__(self):
78   - """
79   - Destructor.
80   - """
81 85 # Return
82 86 return
83 87  
84 88  
85 89 # Private methods
86   - def _parfile(self):
87   - """
88   - Check if parfile exists. If parfile does not exist then create a
89   - default parfile. This kluge avoids shipping the cscript with a parfile.
90   - """
91   - # Set parfile name
92   - parfile = self._name+".par"
93   -
94   - try:
95   - pars = gammalib.GApplicationPars(parfile)
96   - except:
97   - # Signal if parfile was not found
98   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
99   -
100   - # Create default parfile
101   - pars = gammalib.GApplicationPars()
102   - pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Input event list, counts cube, or observation definition XML file"))
103   - pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file"))
104   - pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
105   - pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)"))
106   - pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)"))
107   - pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)"))
108   - pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)"))
109   - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
110   - pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
111   - pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function"))
112   - pars.append(gammalib.GApplicationPar("outfile","f","a","lightcurve.fits","","","Output light curve file"))
113   - pars.append(gammalib.GApplicationPar("tbinalg","s","a","GTI","FILE|LIN|GTI","", "Algorithm for defining time bins"))
114   - pars.append(gammalib.GApplicationPar("tmin","r","a","51544.5","","days","Lightcurve start time [MJD]"))
115   - pars.append(gammalib.GApplicationPar("tmax","r","a","51544.6","","days","Lightcurve stop time [MJD]"))
116   - pars.append(gammalib.GApplicationPar("tbins","i","a","5","","","Number of time bins"))
117   - pars.append(gammalib.GApplicationPar("tbinfile","f","a","tbins.fits","","", "File defining the time binning"))
118   - pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins per light curve bin (0=unbinned)"))
119   - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit of events (TeV)"))
120   - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit of events (TeV)"))
121   - pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
122   - pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method"))
123   - pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)"))
124   - pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)"))
125   - pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
126   - pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
127   - pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
128   - pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value for each bin?"))
129   - pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit for each bin?"))
130   - pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other sky model parameters?"))
131   - pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background model parameters?"))
132   - pars.append_standard()
133   - pars.append(gammalib.GApplicationPar("logfile","f","h","cslightcrv.log","","","Log filename"))
134   - pars.save(parfile)
135   -
136   - # Return
137   - return
138   -
139   - def _create_tbounds(self):
140   -
141   - # Create time bin container
142   - self._tbins = gammalib.GGti()
143   -
144   - # check for binning algorithm
145   - if self["tbinalg"].string() == "FILE":
146   - try:
147   - self._tbins.load(self["tbinfile"].filename())
148   - except:
149   - csv = gammalib.GCsv(self["tbinfile"].filename())
150   - for i in range(csv.nrows()):
151   - tmin = gammalib.GTime()
152   - tmax = gammalib.GTime()
153   - tmin.mjd(csv.real(i,0))
154   - tmax.mjd(csv.real(i,1))
155   - self._tbins.append(tmin,tmax)
156   - elif self["tbinalg"].string() == "LIN":
157   - # Use linear time binning
158   - time_min = self["tmin"].real()
159   - time_max = self["tmax"].real()
160   - nbins = self["tbins"].integer()
161   - step = (time_max - time_min) / float(nbins)
162   - for i in range(nbins):
163   - tmin = gammalib.GTime()
164   - tmin.mjd(time_min + i*step)# ref)
165   - tmax = gammalib.GTime()
166   - tmax.mjd(time_min + (i+1)*step)#ref)
167   - self._tbins.append(tmin,tmax)
168   - elif self["tbinalg"].string() == "GTI":
169   - # Use GTIs of observations
170   - for obs in self._obs:
171   - for i in range(obs.events().gti().size()):
172   - self._tbins.append(obs.events().gti().tstart(i),obs.events().gti().tstop(i))
173   - else:
174   - raise AttributeError("tbinalg=\""+self["tbinalg"].string()+"\" unkown. Must be one of \"FILE\", \"LIN\" or \"GTI\"")
175   -
176   - # Return
177   - return
178   -
179 90 def _get_parameters(self):
180 91 """
181   - Get parameters from parfile and setup the observation.
  92 + Get parameters from parfile.
182 93 """
183 94 # Set observation if not done before
184 95 if self._obs == None or self._obs.size() == 0:
185   - self._require_inobs("cslightcrv::get_parameters()")
  96 + self._require_inobs("cslightcrv::_get_parameters()")
186 97 self._obs = self._get_observations()
187 98  
188 99 # Set models if we have none
... ... @@ -192,46 +103,35 @@ class cslightcrv(ctools.cscript):
192 103 # Get source name
193 104 self._srcname = self["srcname"].string()
194 105  
195   - # Get energy dispersion flag
196   - self._edisp = self["edisp"].boolean()
197   -
198 106 # Get time boundaries
199   - self._create_tbounds()
  107 + self._tbins = self._create_tbounds()
  108 +
  109 + # Set stacked analysis flag if the requested number of energy bins
  110 + # is not zero; otherwise an unbinned analysis will be done.
  111 + self._stacked = False if self["enumbins"].integer() == 0 else True
  112 +
  113 + # Make sure that other parameters are queried now and not later
  114 + self["emin"].real()
  115 + self["emax"].real()
  116 + if self._stacked:
  117 + self["coordsys"].string()
  118 + self["proj"].string()
  119 + self["xref"].real()
  120 + self["yref"].real()
  121 + self["nxpix"].integer()
  122 + self["nypix"].integer()
  123 + self["binsz"].real()
  124 +
  125 + # Do the same for the hidden parameters, just in case
  126 + self["edisp"].boolean()
  127 + self["calc_ulim"].boolean()
  128 + self["calc_ts"].boolean()
  129 + self["fix_bkg"].boolean()
  130 + self["fix_srcs"].boolean()
200 131  
201   - # Unbinned or binned analysis?
202   - self._ebins = self["enumbins"].integer()
203   - if self._ebins == 0:
204   - self._binned = False
205   - else:
206   - self._binned = True
207   -
208   - # Get energy range
209   - self._emin = self["emin"].real()
210   - self._emax = self["emax"].real()
211   -
212   - # Get binning flag
213   - if self._binned:
214   - self._coordsys = self["coordsys"].string()
215   - self._proj = self["proj"].string()
216   - self._xref = self["xref"].real()
217   - self._yref = self["yref"].real()
218   - self._nxpix = self["nxpix"].integer()
219   - self._nypix = self["nypix"].integer()
220   - self._binsz = self["binsz"].real()
221   -
222   - # Read other parameters
223   - self._outfile = self["outfile"].filename()
224   -
225   - # Get other parameeters
226   - self._calc_ulimit = self["calc_ulim"].boolean()
227   - self._calc_ts = self["calc_ts"].boolean()
228   - self._fix_bkg = self["fix_bkg"].boolean()
229   - self._fix_srcs = self["fix_srcs"].boolean()
230   -
231   - # Set some fixed parameters
232   - self._chatter = self["chatter"].integer()
233   - self._clobber = self["clobber"].boolean()
234   - self._debug = self["debug"].boolean()
  132 + # Read ahead output parameters
  133 + if self._read_ahead():
  134 + self["outfile"].filename()
235 135  
236 136 # Write input parameters into logger
237 137 if self._logTerse():
... ... @@ -241,26 +141,101 @@ class cslightcrv(ctools.cscript):
241 141 # Return
242 142 return
243 143  
  144 + def _create_tbounds(self):
  145 + """
  146 + Creates light curve time bins.
244 147  
245   - # Public methods
246   - def run(self):
  148 + The method reads the following user parameters:
  149 + tbinalg: Time binning algorithm
  150 + tbinfile: Time binning file (FITS or ASCII)
  151 + tmin: Start time (MJD)
  152 + tmax: Stop time (MJD)
  153 + tbins: Number of time bins
  154 +
  155 + Returns:
  156 + Light curve bins in form of a GTI.
247 157 """
248   - Run the script.
  158 + # Initialise Good Time Intervals
  159 + gti = gammalib.GGti()
  160 +
  161 + # Get algorithm to use for defining the time intervals
  162 + algorithm = self["tbinalg"].string()
  163 +
  164 + # Handle a FITS or a ASCII file for time bin definition
  165 + if algorithm == "FILE":
  166 +
  167 + # Get the filename
  168 + filename = self["tbinfile"].filename()
  169 +
  170 + # If the file a FITS file then load GTIs
  171 + if filename.is_fits():
  172 + gti.load(filename)
  173 +
  174 + # ... otherwise load file as CSV ASCII file
  175 + csv = gammalib.GCsv(filename)
  176 + for i in range(csv.nrows()):
  177 + tmin = gammalib.GTime()
  178 + tmax = gammalib.GTime()
  179 + tmin.mjd(csv.real(i,0))
  180 + tmax.mjd(csv.real(i,1))
  181 + gti.append(tmin,tmax)
  182 +
  183 + # Handle linear time binning
  184 + elif algorithm == "LIN":
  185 +
  186 + # Get start and stop time and number of time bins
  187 + time_min = self["tmin"].real()
  188 + time_max = self["tmax"].real()
  189 + nbins = self["tbins"].integer()
  190 +
  191 + # Compute time step and setup time intervals
  192 + time_step = (time_max - time_min) / float(nbins)
  193 + for i in range(nbins):
  194 + tmin = gammalib.GTime()
  195 + tmax = gammalib.GTime()
  196 + tmin.mjd(time_min + i *time_step)
  197 + tmax.mjd(time_min + (i+1)*time_step)
  198 + gti.append(tmin,tmax)
  199 +
  200 + # Handle usage of observation GTIs
  201 + elif algorithm == "GTI":
  202 +
  203 + # Append the GTIs of all observations
  204 + for obs in self._obs:
  205 + for i in range(obs.events().gti().size()):
  206 + gti.append(obs.events().gti().tstart(i),
  207 + obs.events().gti().tstop(i))
  208 +
  209 + # ... otherwise raise an exception (this should never occur)
  210 + else:
  211 + raise AttributeError('Paramter tbinalg="'+algorithm+'" unknown. '
  212 + 'Must be one of "FILE", "LIN" or "GTI".')
  213 +
  214 + # Return Good Time Intervals
  215 + return gti
  216 +
  217 + def _get_free_par_names(self):
249 218 """
250   - # Switch screen logging on in debug mode
251   - if self._logDebug():
252   - self._log.cout(True)
  219 + Return list of free parameter names.
  220 +
  221 + Returns:
  222 + List of free parameter names.
  223 + """
  224 + # Initialise list of free parameter names
  225 + names = []
253 226  
254   - # Get parameters
255   - self._get_parameters()
  227 + # Collect list of free parameter names
  228 + for par in self._obs.models()[self._srcname]:
  229 + if par.is_free():
  230 + names.append(par.name())
256 231  
257   - # Write observation into logger
258   - if self._logTerse():
259   - self._log("\n")
260   - self._log.header1("Observation")
261   - self._log(str(self._obs))
262   - self._log("\n")
  232 + # Return names
  233 + return names
263 234  
  235 + def _adjust_model_pars(self):
  236 + """
  237 + Adjust model parameters dependent on user parameters.
  238 + """
264 239 # Write header
265 240 if self._logTerse():
266 241 self._log("\n")
... ... @@ -269,85 +244,341 @@ class cslightcrv(ctools.cscript):
269 244 # Adjust model parameters dependent on input user parameters
270 245 for model in self._obs.models():
271 246  
272   - # Set TS flag for all models to false.
273   - # Source of interest will be set to true later
  247 + # Set TS flag for all models to false. The source of interest
  248 + # will be set to true later
274 249 model.tscalc(False)
275 250  
276 251 # Log model name
277   - if self._logExplicit():
  252 + if self._logNormal():
278 253 self._log.header3(model.name())
279 254  
280 255 # Deal with the source of interest
281 256 if model.name() == self._srcname:
282   - if self._calc_ts:
  257 + if self["calc_ts"].boolean():
283 258 model.tscalc(True)
284 259  
285   - elif self._fix_bkg and not model.classname() == "GModelSky":
  260 + elif (self["fix_bkg"].boolean() and
  261 + model.classname() != "GModelSky"):
286 262 for par in model:
287   - if par.is_free() and self._logExplicit():
288   - self._log(" Fixing \""+par.name()+"\"\n")
289   - par.fix()
290   -
291   - elif self._fix_srcs and model.classname() == "GModelSky":
  263 + if par.is_free():
  264 + par.fix()
  265 + if self._logNormal():
  266 + self._log(gammalib.parformat(par.name()))
  267 + self._log("fixed\n")
  268 +
  269 + elif (self["fix_srcs"].boolean() and
  270 + model.classname() == "GModelSky"):
292 271 for par in model:
293   - if par.is_free() and self._logExplicit():
294   - self._log(" Fixing \""+par.name()+"\"\n")
295   - par.fix()
  272 + if par.is_free():
  273 + par.fix()
  274 + if self._logNormal():
  275 + self._log(gammalib.parformat(par.name()))
  276 + self._log("fixed\n")
296 277  
297   - # Write header
298   - if self._logTerse():
299   - self._log("\n")
300   - self._log.header1("Generate lightcurve")
  278 + # Return
  279 + return
301 280  
302   - # Initialise FITS Table with extension "LIGHTCURVE"
303   - table = gammalib.GFitsBinTable(self._tbins.size())
304   - table.extname("LIGHTCURVE")
  281 + def _create_fits_table(self, results):
  282 + """
  283 + Creates FITS binary table containing light curve results.
  284 +
  285 + Args:
  286 + results: List of result dictionaries.
305 287  
306   - # Add Header for compatibility with gammalib.GMWLSpectrum
307   - table.card("INSTRUME", "CTA", "Name of Instrument")
308   - table.card("TELESCOP", "CTA", "Name of Telescope")
  288 + Returns:
  289 + FITS binary table containing light curve.
  290 + """
  291 + # Determine number of rows in FITS table
  292 + nrows = len(results)
309 293  
310   - # Create FITS table columns
311   - MJD = gammalib.GFitsTableDoubleCol("MJD", self._tbins.size())
312   - MJD.unit("days")
313   - e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self._tbins.size())
314   - e_MJD.unit("days")
  294 + # Create FITS Table with extension "LIGHTCURVE"
  295 + table = gammalib.GFitsBinTable(nrows)
  296 + table.extname("LIGHTCURVE")
315 297  
316   - # Create a FITS column for every free parameter
317   - columns = []
  298 + # Append time columns
  299 + mjd = gammalib.GFitsTableDoubleCol("MJD", nrows)
  300 + e_mjd = gammalib.GFitsTableDoubleCol("e_MJD", nrows)
  301 + mjd.unit("days")
  302 + e_mjd.unit("days")
  303 + for i, result in enumerate(results):
  304 + mjd[i] = result["mjd"]
  305 + e_mjd[i] = result["e_mjd"]
  306 + table.append(mjd)
  307 + table.append(e_mjd)
  308 +
  309 + # Create parameter columns
318 310 for par in self._obs.models()[self._srcname]:
319 311 if par.is_free():
320   - col = gammalib.GFitsTableDoubleCol(par.name(), self._tbins.size())
  312 + name = par.name()
  313 + e_name = "e_"+par.name()
  314 + col = gammalib.GFitsTableDoubleCol(name, nrows)
  315 + e_col = gammalib.GFitsTableDoubleCol(e_name, nrows)
321 316 col.unit(par.unit())
322   - columns.append(col)
323   - e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self._tbins.size())
324 317 e_col.unit(par.unit())
325   - columns.append(e_col)
  318 + for i, result in enumerate(results):
  319 + col[i] = result['values'][name]
  320 + e_col[i] = result['values'][e_name]
  321 + table.append(col)
  322 + table.append(e_col)
  323 +
  324 + # Append Test Statistic column
  325 + ts = gammalib.GFitsTableDoubleCol("TS", nrows)
  326 + for i, result in enumerate(results):
  327 + ts[i] = result["ts"]
  328 + table.append(ts)
  329 +
  330 + # Append upper limit column
  331 + ulimit = gammalib.GFitsTableDoubleCol("UpperLimit", nrows)
  332 + ulimit.unit("ph/cm2/s")
  333 + for i, result in enumerate(results):
  334 + ulimit[i] = result['ulimit']
  335 + table.append(ulimit)
  336 +
  337 + # Return table
  338 + return table
  339 +
  340 + def _compute_ulimit(self, obs):
  341 + """
  342 + Computes upper flux limit.
326 343  
327   - # Create TS and upper limit columns
328   - TSvalues = gammalib.GFitsTableDoubleCol("TS", self._tbins.size())
329   - ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self._tbins.size())
330   - ulim_values.unit("ph/cm2/s")
  344 + Args:
  345 + obs: Observation container.
331 346  
332   - # Loop over energy bins
333   - for i in range(self._tbins.size()):
  347 + Returns:
  348 + Upper flux limit (-1 of not computed).
  349 + """
  350 + # Initialise upper flux limit
  351 + ulimit_value = -1.0
  352 +
  353 + # Perform computation only if requested
  354 + if self["calc_ulim"].boolean():
334 355  
335   - # Log information
336   - if self._logTerse():
337   - self._log("\n")
338   - self._log.header2("Time bin "+str(i))
  356 + # Write header in logger
  357 + if self._logExplicit():
  358 + self._log.header3("Computing upper flux limit")
  359 +
  360 + # Create upper limit object
  361 + ulimit = ctools.ctulimit(obs)
  362 + ulimit["srcname"] = self._srcname
  363 + ulimit["eref"] = 1.0
  364 +
  365 + # Try to run upper limit and catch exceptions
  366 + try:
  367 + ulimit.run()
  368 + ulimit_value = ulimit.flux_ulimit()
  369 + except:
  370 + if self._logTerse():
  371 + self._log("Upper limit flux calculation failed.s\n")
  372 + ulimit_value = -1.0
  373 +
  374 + # Return upper limit
  375 + return ulimit_value
  376 +
  377 + def _bin_observation(self, obs):
  378 + """
  379 + Bin an observation is a binned analysis was requested.
  380 +
  381 + Args:
  382 + obs: Observation container.
  383 +
  384 + Returns:
  385 + Observation container with a binned.
  386 + """
  387 + # Header
  388 + if self._logExplicit():
  389 + self._log.header3("Binning events")
  390 +
  391 + # Bin events
  392 + bin = ctools.ctbin(obs)
  393 + bin["usepnt"] = False
  394 + bin["ebinalg"] = "LOG"
  395 + bin["xref"] = self["xref"].real()
  396 + bin["yref"] = self["yref"].real()
  397 + bin["binsz"] = self["binsz"].real()
  398 + bin["nxpix"] = self["nxpix"].integer()
  399 + bin["nypix"] = self["nypix"].integer()
  400 + bin["enumbins"] = self["enumbins"].integer()
  401 + bin["emin"] = self["emin"].real()
  402 + bin["emax"] = self["emax"].real()
  403 + bin["coordsys"] = self["coordsys"].string()
  404 + bin["proj"] = self["proj"].string()
  405 + bin.run()
  406 +
  407 + # Header
  408 + if self._logExplicit():
  409 + self._log.header3("Creating exposure cube")
  410 +
  411 + # Create exposure cube
  412 + expcube = ctools.ctexpcube(obs)
  413 + expcube["incube"] = "NONE"
  414 + expcube["usepnt"] = False
  415 + expcube["ebinalg"] = "LOG"
  416 + expcube["xref"] = self["xref"].real()
  417 + expcube["yref"] = self["yref"].real()
  418 + expcube["binsz"] = self["binsz"].real()
  419 + expcube["nxpix"] = self["nxpix"].integer()
  420 + expcube["nypix"] = self["nypix"].integer()
  421 + expcube["enumbins"] = self["enumbins"].integer()
  422 + expcube["emin"] = self["emin"].real()
  423 + expcube["emax"] = self["emax"].real()
  424 + expcube["coordsys"] = self["coordsys"].string()
  425 + expcube["proj"] = self["proj"].string()
  426 + expcube.run()
  427 +
  428 + # Header
  429 + if self._logExplicit():
  430 + self._log.header3("Creating point spread function cube")
  431 +
  432 + # Create psf cube
  433 + psfcube = ctools.ctpsfcube(obs)
  434 + psfcube["incube"] = "NONE"
  435 + psfcube["usepnt"] = False
  436 + psfcube["ebinalg"] = "LOG"
  437 + psfcube["xref"] = self["xref"].real()
  438 + psfcube["yref"] = self["yref"].real()
  439 + psfcube["binsz"] = self["binsz"].real()
  440 + psfcube["nxpix"] = self["nxpix"].integer()
  441 + psfcube["nypix"] = self["nypix"].integer()
  442 + psfcube["enumbins"] = self["enumbins"].integer()
  443 + psfcube["emin"] = self["emin"].real()
  444 + psfcube["emax"] = self["emax"].real()
  445 + psfcube["coordsys"] = self["coordsys"].string()
  446 + psfcube["proj"] = self["proj"].string()
  447 + psfcube.run()
  448 +
  449 + # Check if we need to include energy dispersion
  450 + if self["edisp"].boolean():
  451 +
  452 + # Header
  453 + if self._logExplicit():
  454 + self._log.header3("Creating energy dispersion cube")
  455 +
  456 + # Create edisp cube
  457 + edispcube = ctools.ctedispcube(obs)
  458 + edispcube["incube"] = "NONE"
  459 + edispcube["usepnt"] = False
  460 + edispcube["ebinalg"] = "LOG"
  461 + edispcube["xref"] = self["xref"].real()
  462 + edispcube["yref"] = self["yref"].real()
  463 + edispcube["binsz"] = self["binsz"].real()
  464 + edispcube["nxpix"] = self["nxpix"].integer()
  465 + edispcube["nypix"] = self["nypix"].integer()
  466 + edispcube["enumbins"] = self["enumbins"].integer()
  467 + edispcube["emin"] = self["emin"].real()
  468 + edispcube["emax"] = self["emax"].real()
  469 + edispcube["coordsys"] = self["coordsys"].string()
  470 + edispcube["proj"] = self["proj"].string()
  471 + edispcube.run()
  472 +
  473 + # Header
  474 + if self._logExplicit():
  475 + self._log.header3("Creating background cube")
  476 +
  477 + # Create background cube
  478 + bkgcube = ctools.ctbkgcube(obs)
  479 + bkgcube["incube"] = "NONE"
  480 + bkgcube["usepnt"] = False
  481 + bkgcube["ebinalg"] = "LOG"
  482 + bkgcube["xref"] = self["xref"].real()
  483 + bkgcube["yref"] = self["yref"].real()
  484 + bkgcube["binsz"] = self["binsz"].real()
  485 + bkgcube["nxpix"] = self["nxpix"].integer()
  486 + bkgcube["nypix"] = self["nypix"].integer()
  487 + bkgcube["enumbins"] = self["enumbins"].integer()
  488 + bkgcube["emin"] = self["emin"].real()
  489 + bkgcube["emax"] = self["emax"].real()
  490 + bkgcube["coordsys"] = self["coordsys"].string()
  491 + bkgcube["proj"] = self["proj"].string()
  492 + bkgcube.run()
  493 +
  494 + # Retrieve a new oberservation container
  495 + new_obs = bin.obs().copy()
  496 +
  497 + # Get new models
  498 + models = bkgcube.models()
  499 +
  500 + # Set stacked response
  501 + if self["edisp"].boolean():
  502 + new_obs[0].response(expcube.expcube(),
  503 + psfcube.psfcube(),
  504 + edispcube.edispcube(),
  505 + bkgcube.bkgcube())
  506 + else:
  507 + new_obs[0].response(expcube.expcube(),
  508 + psfcube.psfcube(),
  509 + bkgcube.bkgcube())
  510 +
  511 + # Fix background models if required
  512 + if self["fix_bkg"].boolean():
  513 + for model in models:
  514 + if model.classname() != "GModelSky":
  515 + for par in model:
  516 + par.fix()
  517 +
  518 + # Set models for new oberservation container
  519 + new_obs.models(models)
  520 +
  521 + # Return new oberservation container
  522 + return new_obs
  523 +
  524 +
  525 + # Public methods
  526 + def run(self):
  527 + """
  528 + Run the script.
  529 + """
  530 + # Switch screen logging on in debug mode
  531 + if self._logDebug():
  532 + self._log.cout(True)
  533 +
  534 + # Get parameters
  535 + self._get_parameters()
  536 +
  537 + # Write observation into logger
  538 + if self._logTerse():
  539 + self._log("\n")
  540 + self._log.header1("Observation")
  541 + self._log(str(self._obs))
  542 + self._log("\n")
  543 +
  544 + # Adjust model parameters dependent on user parameters
  545 + self._adjust_model_pars()
  546 +
  547 + # Write header
  548 + if self._logTerse():
  549 + self._log("\n")
  550 + self._log.header1("Generate lightcurve")
  551 +
  552 + # Initialise list of result dictionaries
  553 + results = []
  554 +
  555 + # Get source parameters
  556 + pars = self._get_free_par_names()
  557 +
  558 + # Loop over time bins
  559 + for i in range(self._tbins.size()):
339 560  
340 561 # Get time boundaries
341 562 tmin = self._tbins.tstart(i)
342 563 tmax = self._tbins.tstop(i)
343 564  
  565 + # Write time bin into header
  566 + if self._logTerse():
  567 + self._log.header2("MJD "+
  568 + str(tmin.mjd())+"-"+
  569 + str(tmax.mjd()))
  570 +
344 571 # Compute time bin center and time width
345 572 twidth = 0.5 * (tmax - tmin) # in seconds
346 573 tmean = tmin + twidth
347 574  
348   - # Store time as MJD
349   - MJD[i] = tmean.mjd()
350   - e_MJD[i] = twidth / gammalib.sec_in_day # in days
  575 + # Initialise result dictionary
  576 + result = {'mjd': tmean.mjd(),
  577 + 'e_mjd': twidth / gammalib.sec_in_day,
  578 + 'ts': 0.0,
  579 + 'ulimit': 0.0,
  580 + 'pars': pars,
  581 + 'values': {}}
351 582  
352 583 # Log information
353 584 if self._logExplicit():
... ... @@ -355,8 +586,8 @@ class cslightcrv(ctools.cscript):
355 586  
356 587 # Select events
357 588 select = ctools.ctselect(self._obs)
358   - select["emin"] = self._emin
359   - select["emax"] = self._emax
  589 + select["emin"] = self["emin"].real()
  590 + select["emax"] = self["emax"].real()
360 591 select["tmin"] = tmin.convert(select._time_reference())
361 592 select["tmax"] = tmax.convert(select._time_reference())
362 593 select["rad"] = "UNDEFINED"
... ... @@ -367,233 +598,86 @@ class cslightcrv(ctools.cscript):
367 598 # Retrieve observation
368 599 obs = select.obs()
369 600  
370   - # Binned analysis
371   - if self._binned:
372   -
373   - # Header
374   - if self._logTerse():
375   - self._log.header3("Binning events")
376   -
377   - # Bin events
378   - bin = ctools.ctbin(select.obs())
379   - bin["usepnt"] = False
380   - bin["ebinalg"] = "LOG"
381   - bin["xref"] = self._xref
382   - bin["yref"] = self._yref
383   - bin["binsz"] = self._binsz
384   - bin["nxpix"] = self._nxpix
385   - bin["nypix"] = self._nypix
386   - bin["enumbins"] = self._ebins
387   - bin["emin"] = self._emin
388   - bin["emax"] = self._emax
389   - bin["coordsys"] = self._coordsys
390   - bin["proj"] = self._proj
391   - bin.run()
392   -
393   - # Header
394   - if self._logTerse():
395   - self._log.header3("Creating exposure cube")
396   -
397   - # Create exposure cube
398   - expcube = ctools.ctexpcube(select.obs())
399   - expcube["incube"] = "NONE"
400   - expcube["usepnt"] = False
401   - expcube["ebinalg"] = "LOG"
402   - expcube["xref"] = self._xref
403   - expcube["yref"] = self._yref
404   - expcube["binsz"] = self._binsz
405   - expcube["nxpix"] = self._nxpix
406   - expcube["nypix"] = self._nypix
407   - expcube["enumbins"] = self._ebins
408   - expcube["emin"] = self._emin
409   - expcube["emax"] = self._emax
410   - expcube["coordsys"] = self._coordsys
411   - expcube["proj"] = self._proj
412   - expcube.run()
413   -
414   - # Header
415   - if self._logTerse():
416   - self._log.header3("Creating point spread function cube")
417   -
418   - # Create psf cube
419   - psfcube = ctools.ctpsfcube(select.obs())
420   - psfcube["incube"] = "NONE"
421   - psfcube["usepnt"] = False
422   - psfcube["ebinalg"] = "LOG"
423   - psfcube["xref"] = self._xref
424   - psfcube["yref"] = self._yref
425   - psfcube["binsz"] = self._binsz
426   - psfcube["nxpix"] = self._nxpix
427   - psfcube["nypix"] = self._nypix
428   - psfcube["enumbins"] = self._ebins
429   - psfcube["emin"] = self._emin
430   - psfcube["emax"] = self._emax
431   - psfcube["coordsys"] = self._coordsys
432   - psfcube["proj"] = self._proj
433   - psfcube.run()
434   -
435   - # Check if we need to include energy dispersion
436   - if self._edisp:
437   -
438   - # Header
439   - if self._logTerse():
440   - self._log.header3("Creating energy dispersion cube")
441   -
442   - # Create edisp cube
443   - edispcube = ctools.ctedispcube(select.obs())
444   - edispcube["incube"] = "NONE"
445   - edispcube["usepnt"] = False
446   - edispcube["ebinalg"] = "LOG"
447   - edispcube["xref"] = self._xref
448   - edispcube["yref"] = self._yref
449   - edispcube["binsz"] = self._binsz
450   - edispcube["nxpix"] = self._nxpix
451   - edispcube["nypix"] = self._nypix
452   - edispcube["enumbins"] = self._ebins
453   - edispcube["emin"] = self._emin
454   - edispcube["emax"] = self._emax
455   - edispcube["coordsys"] = self._coordsys
456   - edispcube["proj"] = self._proj
457   - edispcube.run()
458   -
459   - # Header
460   - if self._logTerse():
461   - self._log.header3("Creating background cube")
462   -
463   - # Create background cube
464   - bkgcube = ctools.ctbkgcube(select.obs())
465   - bkgcube["incube"] = "NONE"
466   - bkgcube["usepnt"] = False
467   - bkgcube["ebinalg"] = "LOG"
468   - bkgcube["xref"] = self._xref
469   - bkgcube["yref"] = self._yref
470   - bkgcube["binsz"] = self._binsz
471   - bkgcube["nxpix"] = self._nxpix
472   - bkgcube["nypix"] = self._nypix
473   - bkgcube["enumbins"] = self._ebins
474   - bkgcube["emin"] = self._emin
475   - bkgcube["emax"] = self._emax
476   - bkgcube["coordsys"] = self._coordsys
477   - bkgcube["proj"] = self._proj
478   - bkgcube.run()
479   -
480   - # Set new binned observation
481   - obs = bin.obs()
482   -
483   - # Get new models
484   - models = bkgcube.models()
485   -
486   - # Set precomputed binned response
487   - if self._edisp:
488   - obs[0].response(expcube.expcube(), psfcube.psfcube(),
489   - edispcube.edispcube(), bkgcube.bkgcube())
490   - else:
491   - obs[0].response(expcube.expcube(), psfcube.psfcube(),
492   - bkgcube.bkgcube())
493   -
494   - # Fix background models if required
495   - if self._fix_bkg:
496   - for model in models:
497   - if not model.classname() == "GModelSky":
498   - for par in model:
499   - par.fix()
500   -
501   - # Set new models to binned observation
502   - obs.models(models)
  601 + # If a stacked analysis is requested then bin the events
  602 + # and compute the stacked response functions and setup
  603 + # an observation container with a single stacked observation.
  604 + if self._stacked:
  605 + obs = self._bin_observation(obs)
503 606  
504 607 # Header
505   - if self._logTerse():
506   - self._log.header3("Performing fit")
  608 + if self._logExplicit():
  609 + self._log.header3("Fitting the data")
507 610  
508   - # Likelihood
  611 + # Do maximum likelihood model fitting
509 612 like = ctools.ctlike(obs)
510   - like["edisp"] = self._edisp
  613 + like["edisp"] = self["edisp"].boolean()
511 614 like.run()
512 615  
513 616 # Skip bin if no event was present
514 617 if like.obs().logL() == 0.0:
515 618  
516   - # Log information
  619 + # Signal skipping of bin
517 620 if self._logTerse():
518   - self._log("No event in this time bin. Bin is skipped\n")
519   -
520   - # Set all values to 0
521   - for col in columns:
522   - col[i] = 0.0
523   - TSvalues[i] = 0.0
524   - ulim_values[i] = 0.0
525   - continue
  621 + self._log(gammalib.parformat("Warning"))
  622 + self._log("No event in this time bin, skip bin.\n")
526 623  
527   - # Get results
528   - fitted_models = like.obs().models()
529   - source = fitted_models[self._srcname]
  624 + # Set all results to 0
  625 + for par in pars:
  626 + result['values'][par] = 0.0
  627 + result['values']["e_"+par] = 0.0
530 628  
531   - # Calculate Upper Limit
532   - ulimit_value = -1.0
533   - if self._calc_ulimit:
  629 + # Append result
  630 + results.append(result)
534 631  
535   - # Logging information
536   - if self._logTerse():
537   - self._log.header3("Computing upper limit")
538   -
539   - # Create upper limit object
540   - ulimit = ctools.ctulimit(like.obs())
541   - ulimit["srcname"] = self._srcname
542   - ulimit["eref"] = 1.0
543   -
544   - # Try to run upper limit and catch exceptions
545   - try:
546   - ulimit.run()
547   - ulimit_value = ulimit.flux_ulimit()
548   - except:
549   - if self._logTerse():
550   - self._log("Upper limit calculation failed\n")
551   - ulimit_value = -1.0
552   -
553   - # Get TS value
554   - TS = -1.0
555   - if self._calc_ts:
556   - TS = source.ts()
557   -
558   - # Set values for storage
559   - TSvalues[i] = TS
560   -
561   - # Set FITS column values
562   - for col in columns:
563   - if "e_" == col.name()[:2]:
564   - col[i] = source.spectral()[col.name()[2:]].error()
565   - else:
566   - col[i] = source.spectral()[col.name()].value()
567   -
568   - # Store upper limit value if available
569   - if ulimit_value > 0.0:
570   - ulim_values[i] = ulimit_value
  632 + # Continue with next time bin
  633 + continue
571 634  
572   - # Log information
573   - if self._logExplicit():
574   - self._log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd()))
575   - for col in columns:
576   - if "e_" == col.name()[:2]:
577   - continue
578   - value = source.spectral()[col.name()].value()
579   - error = source.spectral()[col.name()].error()
580   - unit = source.spectral()[col.name()].unit()
581   - self._log(" > "+col.name()+": "+str(value)+" +- "+str(error)+" "+unit+"\n")
582   - if self._calc_ts and TSvalues[i] > 0.0:
583   - self._log(" > TS = "+str(TS)+" \n")
584   - if self._calc_ulimit and ulim_values[i] > 0.0:
585   - self._log(" > UL = "+str(ulim_values[i])+" [ph/cm2/s]")
586   - self._log("\n")
  635 + # Retrieve model fitting results for source of interest
  636 + source = like.obs().models()[self._srcname]
587 637  
588   - # Append filles columns to fits table
589   - table.append(MJD)
590   - table.append(e_MJD)
591   - for col in columns:
592   - table.append(col)
593   - table.append(TSvalues)
594   - table.append(ulim_values)
  638 + # Extract parameter values
  639 + for par in pars:
  640 + result['values'][par] = source.spectral()[par].value()
  641 + result['values']["e_"+par] = source.spectral()[par].error()
595 642  
596   - # Create the FITS file now
  643 + # Calculate upper limit (-1 if not computed)
  644 + ulimit_value = self._compute_ulimit(like.obs())
  645 + if ulimit_value > 0.0:
  646 + result['ulimit'] = ulimit_value
  647 +
  648 + # Extract Test Statistic value
  649 + if self["calc_ts"].boolean():
  650 + result['ts'] = source.ts()
  651 +
  652 + # Append result to list of dictionaries
  653 + results.append(result)
  654 +
  655 + # Log results for this time bin
  656 + if self._logNormal():
  657 + self._log.header3("Results")
  658 + pars = self._get_free_par_names()
  659 + for par in pars:
  660 + value = source.spectral()[par].value()
  661 + error = source.spectral()[par].error()
  662 + unit = source.spectral()[par].unit()
  663 + self._log(gammalib.parformat(par))
  664 + self._log(str(value))
  665 + self._log(" +/- ")
  666 + self._log(str(error))
  667 + self._log(" ")
  668 + self._log(unit)
  669 + self._log("\n")
  670 + if result['ulimit'] > 0.0:
  671 + self._log(gammalib.parformat("Upper flux limit"))
  672 + self._log(str(result['ulimit'])+" ph/cm2/s\n")
  673 + if self["calc_ts"].boolean():
  674 + self._log(gammalib.parformat("Test Statistic"))
  675 + self._log(str(result['ts'])+"\n")
  676 +
  677 + # Create FITS table from results
  678 + table = self._create_fits_table(results)
  679 +
  680 + # Create FITS file and append FITS table to FITS file
597 681 self._fits = gammalib.GFits()
598 682 self._fits.append(table)
599 683  
... ... @@ -604,6 +688,9 @@ class cslightcrv(ctools.cscript):
604 688 """
605 689 Execute the script.
606 690 """
  691 + # Open logfile
  692 + self._logFileOpen()
  693 +
607 694 # Read ahead output parameters
608 695 self._read_ahead(True)
609 696  
... ... @@ -618,35 +705,37 @@ class cslightcrv(ctools.cscript):
618 705  
619 706 def save(self):
620 707 """
621   - Save lightcurve.
  708 + Save light curve.
622 709 """
623 710 # Write header
624 711 if self._logTerse():
625   - self._log("\n")
626   - self._log.header1("Save lightcurve")
  712 + self._log('\n')
  713 + self._log.header1('Save light curve')
627 714  
628   - # Get outmap parameter
629   - self._outfile = self["outfile"].filename()
  715 + # Get light curve filename
  716 + outfile = self["outfile"].filename()
630 717  
631 718 # Continue only filename and residual map are valid
632   - if self._outfile != "NONE" and self._fits != None:
  719 + if self._fits != None:
633 720  
634 721 # Log file name
635 722 if self._logTerse():
636   - self._log("Save lightcurve into \""+self._outfile+"\".\n")
  723 + self._log(gammalib.parformat("Light curve file"))
  724 + self._log(outfile.url())
  725 + self._log("\n")
637 726  
638 727 # Save spectrum
639   - self._fits.saveto(self._outfile, self._clobber)
  728 + self._fits.saveto(outfile, self._clobber())
640 729  
641 730 # Return
642 731 return
643 732  
644 733 def lightcurve(self):
645 734 """
646   - Return lightcurve FITS file.
  735 + Return light curve FITS file.
647 736  
648 737 Returns:
649   - FITS file containing lightcurve.
  738 + FITS file containing light curve.
650 739 """
651 740 # Return
652 741 return self._fits
... ... @@ -670,8 +759,5 @@ if __name__ == '__main__':
670 759 # Create instance of application
671 760 app = cslightcrv(sys.argv)
672 761  
673   - # Open logfile
674   - app._logFileOpen()
675   -
676 762 # Execute application
677 763 app.execute()
... ...
cscripts/csmodelinfo.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: csmodelinfo.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the csmodelinfo script #
15 15 # #
... ...
cscripts/csmodelinfo.py
  View file @31f23bc
... ... @@ -36,7 +36,7 @@ class csmodelinfo(ctools.cscript):
36 36 with DS9.
37 37 """
38 38  
39   - # Constructors and destructors
  39 + # Constructor
40 40 def __init__(self, *argv):
41 41 """
42 42 Constructor.
... ... @@ -47,7 +47,7 @@ class csmodelinfo(ctools.cscript):
47 47  
48 48 # Initialise class members
49 49 self._models = gammalib.GModels()
50   - self._ds9file = "NONE"
  50 + self._ds9file = gammalib.GFilename("NONE")
51 51 self._pnt_type = ""
52 52 self._pnt_mark_size = 12
53 53 self._show_labels = True
... ... @@ -59,27 +59,11 @@ class csmodelinfo(ctools.cscript):
59 59 self._show_ext_type = True
60 60 self._free_color = "green"
61 61 self._fixed_color = "magenta"
62   - self._read_ahead = False
63 62  
64   - # Initialise application
65   - if len(argv) == 0:
66   - ctools.cscript.__init__(self, self._name, self._version)
67   - elif len(argv) == 1:
68   - ctools.cscript.__init__(self, self._name, self._version, *argv)
69   - else:
70   - raise TypeError("Invalid number of arguments given.")
  63 + # Initialise application by calling the appropriate class
  64 + # constructor.
  65 + self._init_cscript(argv)
71 66  
72   - # Set logger properties
73   - self._log_header()
74   - self._log.date(True)
75   -
76   - # Return
77   - return
78   -
79   - def __del__(self):
80   - """
81   - Destructor.
82   - """
83 67 # Return
84 68 return
85 69  
... ... @@ -106,7 +90,7 @@ class csmodelinfo(ctools.cscript):
106 90 self._fixed_color = self["fixed_color"].string()
107 91  
108 92 # Read ahead DS9 filename
109   - if self._read_ahead:
  93 + if self._read_ahead():
110 94 self._ds9file = self["ds9file"].filename()
111 95  
112 96 # Write input parameters into logger
... ... @@ -500,8 +484,11 @@ class csmodelinfo(ctools.cscript):
500 484 """
501 485 Execute the script.
502 486 """
503   - # Set read ahead flag
504   - self._read_ahead = True
  487 + # Open logfile
  488 + self._logFileOpen()
  489 +
  490 + # Read ahead output parameters
  491 + self._read_ahead(True)
505 492  
506 493 # Run the script
507 494 self.run()
... ... @@ -521,8 +508,5 @@ if __name__ == '__main__':
521 508 # Create instance of application
522 509 app = csmodelinfo(sys.argv)
523 510  
524   - # Open logfile
525   - app._logFileOpen()
526   -
527 511 # Execute application
528 512 app.execute()
... ...
cscripts/csmodelmerge.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: csmodelmerge.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the csmodelmerge script #
15 15 # #
... ...
cscripts/csmodelmerge.py
  View file @31f23bc
... ... @@ -35,7 +35,7 @@ class csmodelmerge(ctools.cscript):
35 35 a single model definition XML file.
36 36 """
37 37  
38   - # Constructors and destructors
  38 + # Constructor
39 39 def __init__(self, *argv):
40 40 """
41 41 Constructor.
... ... @@ -47,28 +47,11 @@ class csmodelmerge(ctools.cscript):
47 47 # Initialise class members
48 48 self._files = None
49 49 self._models = gammalib.GModels()
50   - self._outmodel = "NONE"
51   - self._read_ahead = False
52   -
53   - # Initialise application
54   - if len(argv) == 0:
55   - ctools.cscript.__init__(self, self._name, self._version)
56   - elif len(argv) == 1:
57   - ctools.cscript.__init__(self, self._name, self._version, *argv)
58   - else:
59   - raise TypeError("Invalid number of arguments given.")
60   -
61   - # Set logger properties
62   - self._log_header()
63   - self._log.date(True)
64 50  
65   - # Return
66   - return
  51 + # Initialise application by calling the appropriate class
  52 + # constructor.
  53 + self._init_cscript(argv)
67 54  
68   - def __del__(self):
69   - """
70   - Destructor.
71   - """
72 55 # Return
73 56 return
74 57  
... ... @@ -105,8 +88,8 @@ class csmodelmerge(ctools.cscript):
105 88 raise RuntimeError(msg)
106 89  
107 90 # Read ahead output filename
108   - if self._read_ahead:
109   - self.outmodel = self["outmodel"].filename()
  91 + if self._read_ahead():
  92 + self["outmodel"].filename()
110 93  
111 94 # Write input parameters into logger
112 95 if self._logTerse():
... ... @@ -178,16 +161,16 @@ class csmodelmerge(ctools.cscript):
178 161 self._log.header1("Save models")
179 162  
180 163 # Get output filename in case it was not read ahead
181   - self._outmodel = self["outmodel"].filename()
  164 + outmodel = self["outmodel"].filename()
182 165  
183 166 # Log filename
184 167 if self._logTerse():
185 168 self._log(gammalib.parformat("Model definition XML file"))
186   - self._log(self._outmodel.url())
  169 + self._log(outmodel.url())
187 170 self._log("\n")
188 171  
189 172 # Save models
190   - self._models.save(self._outmodel)
  173 + self._models.save(outmodel)
191 174  
192 175 # Return
193 176 return
... ... @@ -196,8 +179,11 @@ class csmodelmerge(ctools.cscript):
196 179 """
197 180 Execute the script.
198 181 """
199   - # Set read ahead flag
200   - self._read_ahead = True
  182 + # Open logfile
  183 + self._logFileOpen()
  184 +
  185 + # Read ahead output parameters
  186 + self._read_ahead(True)
201 187  
202 188 # Run the script
203 189 self.run()
... ... @@ -217,8 +203,5 @@ if __name__ == '__main__':
217 203 # Create instance of application
218 204 app = csmodelmerge(sys.argv)
219 205  
220   - # Open logfile
221   - app._logFileOpen()
222   -
223 206 # Execute application
224 207 app.execute()
... ...
cscripts/csobs2caldb.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: csobs2caldb.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the csobs2caldb script #
15 15 # #
... ...
cscripts/csobs2caldb.py
  View file @31f23bc
... ... @@ -40,7 +40,7 @@ class csobs2caldb(ctools.cscript):
40 40 index of any observation using the hidden "index" parameter.
41 41 """
42 42  
43   - # Constructors and destructors
  43 + # Constructor
44 44 def __init__(self, *argv):
45 45 """
46 46 Constructor.
... ... @@ -59,28 +59,13 @@ class csobs2caldb(ctools.cscript):
59 59 self._base_dir = ""
60 60 self._cal_dir = ""
61 61  
62   - # Initialise application
63   - if len(argv) == 0:
64   - ctools.cscript.__init__(self, self._name, self._version)
65   - elif len(argv) == 1:
66   - ctools.cscript.__init__(self, self._name, self._version, *argv)
67   - else:
68   - raise TypeError("Invalid number of arguments given.")
69   -
70   - # Set logger properties
71   - self._log_header()
72   - self._log.date(True)
  62 + # Initialise application by calling the appropriate class
  63 + # constructor.
  64 + self._init_cscript(argv)
73 65  
74 66 # Return
75 67 return
76 68  
77   - def __del__(self):
78   - """
79   - Destructor.
80   - """
81   - # Return
82   - return
83   -
84 69  
85 70 # Private methods
86 71 def _get_parameters(self):
... ... @@ -395,6 +380,9 @@ class csobs2caldb(ctools.cscript):
395 380 """
396 381 Execute the script.
397 382 """
  383 + # Open logfile
  384 + self._logFileOpen()
  385 +
398 386 # Run the script
399 387 self.run()
400 388  
... ... @@ -413,8 +401,5 @@ if __name__ == '__main__':
413 401 # Create instance of application
414 402 app = csobs2caldb(sys.argv)
415 403  
416   - # Open logfile
417   - app._logFileOpen()
418   -
419 404 # Execute application
420 405 app.execute()
... ...
cscripts/csobsdef.py
  View file @31f23bc
... ... @@ -72,7 +72,7 @@ class csobsdef(ctools.cscript):
72 72 IRF in the "dummy" calibration database.
73 73 """
74 74  
75   - # Constructors and destructors
  75 + # Constructor
76 76 def __init__(self, *argv):
77 77 """
78 78 Constructor.
... ... @@ -85,31 +85,14 @@ class csobsdef(ctools.cscript):
85 85 self._obs = gammalib.GObservations()
86 86 self._inpnt = ""
87 87 self._tmin = 0.0
88   - self._outobs = "NONE"
89 88 self._chatter = 2
90 89 self._clobber = True
91 90 self._debug = False
92   - self._read_ahead = False
93 91  
94   - # Initialise application
95   - if len(argv) == 0:
96   - ctools.cscript.__init__(self, self._name, self._version)
97   - elif len(argv) ==1:
98   - ctools.cscript.__init__(self, self._name, self._version, *argv)
99   - else:
100   - raise TypeError("Invalid number of arguments given.")
  92 + # Initialise application by calling the appropriate class
  93 + # constructor.
  94 + self._init_cscript(argv)
101 95  
102   - # Set logger properties
103   - self._log_header()
104   - self._log.date(True)
105   -
106   - # Return
107   - return
108   -
109   - def __del__(self):
110   - """
111   - Destructor.
112   - """
113 96 # Return
114 97 return
115 98  
... ... @@ -123,8 +106,8 @@ class csobsdef(ctools.cscript):
123 106 self._inpnt = self["inpnt"].filename()
124 107  
125 108 # Read ahead parameters
126   - if self._read_ahead:
127   - self._outobs = self["outobs"].filename()
  109 + if self._read_ahead():
  110 + self["outobs"].filename()
128 111  
129 112 # Set some fixed parameters
130 113 self._chatter = self["chatter"].integer()
... ... @@ -347,19 +330,19 @@ class csobsdef(ctools.cscript):
347 330 self._log.header1("Save observation definition XML file")
348 331  
349 332 # Get output filename in case it was not read ahead
350   - self._outobs = self["outobs"].filename()
  333 + outobs = self["outobs"].filename()
351 334  
352 335 # Check if observation definition XML file is valid
353   - if self._outobs != "NONE":
  336 + if outobs.url() != "NONE":
354 337  
355 338 # Log filename
356 339 if self._logTerse():
357 340 self._log(gammalib.parformat("Observation XML file"))
358   - self._log(self._outobs.url())
  341 + self._log(outobs.url())
359 342 self._log("\n")
360 343  
361 344 # Save observation definition XML file
362   - self._obs.save(self._outobs)
  345 + self._obs.save(outobs)
363 346  
364 347 # Return
365 348 return
... ... @@ -368,8 +351,11 @@ class csobsdef(ctools.cscript):
368 351 """
369 352 Execute the script.
370 353 """
371   - # Set read ahead flag
372   - self._read_ahead = True
  354 + # Open logfile
  355 + self._logFileOpen()
  356 +
  357 + # Read ahead output parameters
  358 + self._read_ahead(True)
373 359  
374 360 # Run the script
375 361 self.run()
... ... @@ -389,8 +375,5 @@ if __name__ == '__main__':
389 375 # Create instance of application
390 376 app = csobsdef(sys.argv)
391 377  
392   - # Open logfile
393   - app._logFileOpen()
394   -
395 378 # Execute application
396 379 app.execute()
... ...
cscripts/csobsinfo.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: csobsinfo.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the csobsinfo script #
15 15 # #
... ...
cscripts/csobsinfo.py
  View file @31f23bc
... ... @@ -31,7 +31,7 @@ class csobsinfo(ctools.cscript):
31 31 Shows the content of an observation container.
32 32 """
33 33  
34   - # Constructors and destructors
  34 + # Constructor
35 35 def __init__(self, *argv):
36 36 """
37 37 Constructor.
... ... @@ -43,7 +43,6 @@ class csobsinfo(ctools.cscript):
43 43 # Initialise class members
44 44 self._obj_dir = None
45 45 self._compute_offset = False
46   - self._ds9file = ""
47 46 self._offsets = []
48 47 self._zeniths = []
49 48 self._azimuths = []
... ... @@ -51,36 +50,17 @@ class csobsinfo(ctools.cscript):
51 50 self._pnt_dec = []
52 51 self._ebounds = gammalib.GEbounds()
53 52 self._gti = gammalib.GGti()
54   - self._read_ahead = False
55   - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
56   - self._obs = argv[0]
57   - argv = argv[1:]
58   - else:
59   - self._obs = gammalib.GObservations()
60   - self._obs.clear()
61 53  
62   - # Initialise application
63   - if len(argv) == 0:
64   - ctools.cscript.__init__(self, self._name, self._version)
65   - elif len(argv) == 1:
66   - ctools.cscript.__init__(self, self._name, self._version, *argv)
67   - else:
68   - raise TypeError("Invalid number of arguments given.")
  54 + # Initialise observation container from constructor arguments.
  55 + self._obs = self._set_input_obs(argv)
69 56  
70   - # Set logger properties
71   - self._log_header()
72   - self._log.date(True)
  57 + # Initialise application by calling the appropriate class
  58 + # constructor.
  59 + self._init_cscript(argv)
73 60  
74 61 # Return
75 62 return
76 63  
77   - def __del__(self):
78   - """
79   - Destructor.
80   - """
81   - # Return
82   - return
83   -
84 64  
85 65 # Private methods
86 66 def _get_parameters(self):
... ... @@ -103,8 +83,8 @@ class csobsinfo(ctools.cscript):
103 83 self._obj_dir.radec_deg(ra,dec)
104 84  
105 85 # Read ahead DS9 filename
106   - if self._read_ahead:
107   - self._ds9file = self["ds9file"].filename()
  86 + if self._read_ahead():
  87 + self["ds9file"].filename()
108 88  
109 89 # Write input parameters into logger
110 90 if self._logTerse():
... ... @@ -387,19 +367,19 @@ class csobsinfo(ctools.cscript):
387 367 self._log.header1("Save pointings in DS9 file")
388 368  
389 369 # Get output filename in case it was not read ahead
390   - self._ds9file = self["ds9file"].filename()
  370 + ds9file = self["ds9file"].filename()
391 371  
392 372 # Check if DS9 file is valid
393   - if self._ds9file != "NONE":
  373 + if ds9file.url() != "NONE":
394 374  
395 375 # Log filename
396 376 if self._logTerse():
397 377 self._log(gammalib.parformat("DS9 filename"))
398   - self._log(self._ds9file.url())
  378 + self._log(ds9file.url())
399 379 self._log("\n")
400 380  
401 381 # Open file
402   - f = open(self._ds9file.url(),"w")
  382 + f = open(ds9file.url(),"w")
403 383  
404 384 # Write coordinate system
405 385 f.write("fk5\n")
... ... @@ -455,8 +435,11 @@ class csobsinfo(ctools.cscript):
455 435 """
456 436 Execute the script.
457 437 """
458   - # Set read ahead flag
459   - self._read_ahead = True
  438 + # Open logfile
  439 + self._logFileOpen()
  440 +
  441 + # Read ahead output parameters
  442 + self._read_ahead(True)
460 443  
461 444 # Run the script
462 445 self.run()
... ... @@ -476,8 +459,5 @@ if __name__ == '__main__':
476 459 # Create instance of application
477 460 app = csobsinfo(sys.argv)
478 461  
479   - # Open logfile
480   - app._logFileOpen()
481   -
482 462 # Execute application
483 463 app.execute()
... ...
cscripts/cspull.py
  View file @31f23bc
... ... @@ -38,7 +38,7 @@ class cspull(ctools.cscript):
38 38 script behaves just as a regular ctool.
39 39 """
40 40  
41   - # Constructors and destructors
  41 + # Constructor
42 42 def __init__(self, *argv):
43 43 """
44 44 Constructor.
... ... @@ -48,7 +48,6 @@ class cspull(ctools.cscript):
48 48 self._version = "1.1.0"
49 49  
50 50 # Initialise some members
51   - self._obs = None
52 51 self._model = None
53 52 self._inmodel = None
54 53 self._edisp = False
... ... @@ -62,25 +61,13 @@ class cspull(ctools.cscript):
62 61 self._proj = "TAN"
63 62 self._log_clients = False
64 63  
65   - # Initialise application
66   - if len(argv) == 0:
67   - ctools.cscript.__init__(self, self._name, self._version)
68   - elif len(argv) ==1:
69   - ctools.cscript.__init__(self, self._name, self._version, *argv)
70   - else:
71   - raise TypeError("Invalid number of arguments given.")
72   -
73   - # Set logger properties
74   - self._log_header()
75   - self._log.date(True)
  64 + # Initialise observation container from constructor arguments.
  65 + self._obs = self._set_input_obs(argv)
76 66  
77   - # Return
78   - return
  67 + # Initialise application by calling the appropriate class
  68 + # constructor.
  69 + self._init_cscript(argv)
79 70  
80   - def __del__(self):
81   - """
82   - Destructor.
83   - """
84 71 # Return
85 72 return
86 73  
... ... @@ -310,10 +297,10 @@ class cspull(ctools.cscript):
310 297  
311 298 # Fit model
312 299 if self._profile:
313   - models = _obs.models()
  300 + models = self._obs.models()
314 301 for i in range(models.size()):
315 302 model_name = models[i].name()
316   - like = obsutils.cterror(_obs, model_name,
  303 + like = obsutils.cterror(self._obs, model_name,
317 304 log=self._log_clients,
318 305 debug=self._debug,
319 306 chatter=self._chatter)
... ... @@ -442,6 +429,9 @@ class cspull(ctools.cscript):
442 429 """
443 430 Execute the script.
444 431 """
  432 + # Open logfile
  433 + self._logFileOpen()
  434 +
445 435 # Run the script
446 436 self.run()
447 437  
... ... @@ -467,8 +457,5 @@ if __name__ == '__main__':
467 457 # Create instance of application
468 458 app = cspull(sys.argv)
469 459  
470   - # Open logfile
471   - app._logFileOpen()
472   -
473 460 # Execute application
474 461 app.execute()
... ...
cscripts/csresmap.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: csresmap.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the csresmap script #
15 15 # #
... ...
cscripts/csresmap.py
  View file @31f23bc
... ... @@ -36,7 +36,7 @@ class csresmap(ctools.cscript):
36 36 behaves just as a regular ctool.
37 37 """
38 38  
39   - # Constructors and destructors
  39 + # Constructor
40 40 def __init__(self, *argv):
41 41 """
42 42 Constructor.
... ... @@ -49,7 +49,6 @@ class csresmap(ctools.cscript):
49 49 self._algorithm = "SUB"
50 50 self._resmap = None
51 51 self._modcube = "NONE"
52   - self._outmap = "NONE"
53 52 self._use_maps = False
54 53 self._skip_binning = False
55 54 self._publish = False
... ... @@ -74,35 +73,16 @@ class csresmap(ctools.cscript):
74 73 self._debug = False
75 74 self._log_clients = False
76 75  
77   - # Initialise observation
78   - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
79   - self._obs = argv[0]
80   - argv = argv[1:]
81   - else:
82   - self._obs = gammalib.GObservations()
83   -
84   - # Initialise application
85   - if len(argv) == 0:
86   - ctools.cscript.__init__(self, self._name, self._version)
87   - elif len(argv) ==1:
88   - ctools.cscript.__init__(self, self._name, self._version, *argv)
89   - else:
90   - raise TypeError("Invalid number of arguments given.")
  76 + # Initialise observation container from constructor arguments.
  77 + self._obs = self._set_input_obs(argv)
91 78  
92   - # Set logger properties
93   - self._log_header()
94   - self._log.date(True)
  79 + # Initialise application by calling the appropriate class
  80 + # constructor.
  81 + self._init_cscript(argv)
95 82  
96 83 # Return
97 84 return
98 85  
99   - def __del__(self):
100   - """
101   - Destructor.
102   - """
103   - # Return
104   - return
105   -
106 86  
107 87 # Private methods
108 88 def _get_parameters(self):
... ... @@ -177,7 +157,7 @@ class csresmap(ctools.cscript):
177 157  
178 158 # Read ahead output parameters
179 159 if (self._read_ahead()):
180   - self._outmap = self["outmap"].filename()
  160 + self["outmap"].filename()
181 161  
182 162 # Write input parameters into logger
183 163 if self._logTerse():
... ... @@ -318,6 +298,9 @@ class csresmap(ctools.cscript):
318 298 """
319 299 Execute the script.
320 300 """
  301 + # Open logfile
  302 + self._logFileOpen()
  303 +
321 304 # Read ahead output parameters
322 305 self._read_ahead(True)
323 306  
... ... @@ -340,17 +323,19 @@ class csresmap(ctools.cscript):
340 323 self._log.header1("Save residual map")
341 324  
342 325 # Get outmap parameter
343   - self._outmap = self["outmap"].filename()
  326 + outmap = self["outmap"].filename()
344 327  
345 328 # Continue only filename and residual map are valid
346   - if self._outmap != "NONE" and self._resmap != None:
  329 + if self._resmap != None:
347 330  
348 331 # Log file name
349 332 if self._logTerse():
350   - self._log("Save residual map into \""+self._outmap+"\".\n")
  333 + self._log(gammalib.parformat("Residual map file"))
  334 + self._log(outmap.url())
  335 + self._log("\n")
351 336  
352 337 # Save residual map
353   - self._resmap.save(self._outmap, self._clobber)
  338 + self._resmap.save(outmap, self._clobber)
354 339  
355 340 # Return
356 341 return
... ... @@ -405,8 +390,5 @@ if __name__ == '__main__':
405 390 # Create instance of application
406 391 app = csresmap(sys.argv)
407 392  
408   - # Open logfile
409   - app._logFileOpen()
410   -
411 393 # Execute application
412 394 app.execute()
... ...
cscripts/cssens.py
  View file @31f23bc
... ... @@ -41,7 +41,7 @@ class cssens(ctools.cscript):
41 41 the likelihood difference between fitting with and without the test source.
42 42 """
43 43  
44   - # Constructors and destructors
  44 + # Constructor
45 45 def __init__(self, *argv):
46 46 """
47 47 Constructor.
... ... @@ -54,7 +54,7 @@ class cssens(ctools.cscript):
54 54 self._obs = gammalib.GObservations()
55 55 self._ebounds = gammalib.GEbounds()
56 56 self._srcname = ""
57   - self._outfile = ""
  57 + self._outfile = gammalib.GFilename()
58 58 self._ra = None
59 59 self._dec = None
60 60 self._edisp = False
... ... @@ -71,28 +71,13 @@ class cssens(ctools.cscript):
71 71 self._log_clients = False
72 72 self._debug = False
73 73  
74   - # Initialise application
75   - if len(argv) == 0:
76   - ctools.cscript.__init__(self, self._name, self._version)
77   - elif len(argv) ==1:
78   - ctools.cscript.__init__(self, self._name, self._version, *argv)
79   - else:
80   - raise TypeError("Invalid number of arguments given.")
81   -
82   - # Set logger properties
83   - self._log_header()
84   - self._log.date(True)
  74 + # Initialise application by calling the appropriate class
  75 + # constructor.
  76 + self._init_cscript(argv)
85 77  
86 78 # Return
87 79 return
88 80  
89   - def __del__(self):
90   - """
91   - Destructor.
92   - """
93   - # Return
94   - return
95   -
96 81  
97 82 # Private methods
98 83 def _get_parameters(self):
... ... @@ -734,6 +719,9 @@ class cssens(ctools.cscript):
734 719 """
735 720 Execute the script.
736 721 """
  722 + # Open logfile
  723 + self._logFileOpen()
  724 +
737 725 # Run the script
738 726 self.run()
739 727  
... ... @@ -750,8 +738,5 @@ if __name__ == '__main__':
750 738 # Create instance of application
751 739 app = cssens(sys.argv)
752 740  
753   - # Open logfile
754   - app._logFileOpen()
755   -
756 741 # Run application
757 742 app.execute()
... ...
cscripts/csspec.py
  View file @31f23bc
... ... @@ -51,7 +51,6 @@ class csspec(ctools.cscript):
51 51 self._binned_mode = False
52 52 self._srcname = ""
53 53 self._edisp = False
54   - self._outfile = "NONE"
55 54 self._calc_ulimit = True
56 55 self._calc_ts = True
57 56 self._fix_bkg = False
... ... @@ -60,24 +59,12 @@ class csspec(ctools.cscript):
60 59 self._clobber = True
61 60 self._debug = False
62 61  
63   - # Initialise observation
64   - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
65   - self._obs = argv[0]
66   - argv = argv[1:]
67   - else:
68   - self._obs = gammalib.GObservations()
69   -
70   - # Initialise application
71   - if len(argv) == 0:
72   - ctools.cscript.__init__(self, self._name, self._version)
73   - elif len(argv) ==1:
74   - ctools.cscript.__init__(self, self._name, self._version, *argv)
75   - else:
76   - raise TypeError("Invalid number of arguments given.")
  62 + # Initialise observation container from constructor arguments.
  63 + self._obs = self._set_input_obs(argv)
77 64  
78   - # Set logger properties
79   - self._log_header()
80   - self._log.date(True)
  65 + # Initialise application by calling the appropriate class
  66 + # constructor.
  67 + self._init_cscript(argv)
81 68  
82 69 # Return
83 70 return
... ... @@ -132,8 +119,8 @@ class csspec(ctools.cscript):
132 119 self._debug = self["debug"].boolean()
133 120  
134 121 # Read ahead output parameters
135   - if (self._read_ahead()):
136   - self._outfile = self["outfile"].filename()
  122 + if self._read_ahead():
  123 + self["outfile"].filename()
137 124  
138 125 # Write input parameters into logger
139 126 if self._logTerse():
... ... @@ -494,6 +481,9 @@ class csspec(ctools.cscript):
494 481 """
495 482 Execute the script.
496 483 """
  484 + # Open logfile
  485 + self._logFileOpen()
  486 +
497 487 # Read ahead output parameters
498 488 self._read_ahead(True)
499 489  
... ... @@ -516,17 +506,19 @@ class csspec(ctools.cscript):
516 506 self._log.header1("Save spectrum")
517 507  
518 508 # Get outmap parameter
519   - self._outfile = self["outfile"].filename()
  509 + outfile = self["outfile"].filename()
520 510  
521 511 # Continue only filename and residual map are valid
522   - if self._outfile != "NONE" and self._fits != None:
  512 + if self._fits != None:
523 513  
524 514 # Log file name
525 515 if self._logTerse():
526   - self._log("Save spectrum into \""+self._outfile+"\".\n")
  516 + self._log(gammalib.parformat("Spectrum file"))
  517 + self._log(outfile.url())
  518 + self._log("\n")
527 519  
528 520 # Save spectrum
529   - self._fits.saveto(self._outfile, self._clobber)
  521 + self._fits.saveto(outfile, self._clobber)
530 522  
531 523 # Return
532 524 return
... ... @@ -560,8 +552,5 @@ if __name__ == '__main__':
560 552 # Create instance of application
561 553 app = csspec(sys.argv)
562 554  
563   - # Open logfile
564   - app._logFileOpen()
565   -
566 555 # Execute application
567 556 app.execute()
... ...
cscripts/cstsdist.py
  View file @31f23bc
... ... @@ -39,7 +39,7 @@ class cstsdist(ctools.cscript):
39 39 script behaves just as a regular ctool.
40 40 """
41 41  
42   - # Constructors and destructors
  42 + # Constructor
43 43 def __init__(self, *argv):
44 44 """
45 45 Constructor.
... ... @@ -55,34 +55,19 @@ class cstsdist(ctools.cscript):
55 55 self._enumbins = 0
56 56 self._npix = 0
57 57 self._binsz = 0.0
58   - self._outfile = "ts.dat"
  58 + self._outfile = gammalib.GFilename("ts.dat")
59 59 self._ntrials = 10
60 60 self._edisp = False
61 61 self._debug = False
62 62 self._log_clients = False
63 63  
64   - # Initialise application
65   - if len(argv) == 0:
66   - ctools.cscript.__init__(self, self._name, self._version)
67   - elif len(argv) ==1:
68   - ctools.cscript.__init__(self, self._name, self._version, *argv)
69   - else:
70   - raise TypeError("Invalid number of arguments given.")
71   -
72   - # Set logger properties
73   - self._log_header()
74   - self._log.date(True)
  64 + # Initialise application by calling the appropriate class
  65 + # constructor.
  66 + self._init_cscript(argv)
75 67  
76 68 # Return
77 69 return
78 70  
79   - def __del__(self):
80   - """
81   - Destructor.
82   - """
83   - # Return
84   - return
85   -
86 71  
87 72 # Private methods
88 73 def _get_parameters(self):
... ... @@ -463,6 +448,9 @@ class cstsdist(ctools.cscript):
463 448 """
464 449 Execute the script.
465 450 """
  451 + # Open logfile
  452 + self._logFileOpen()
  453 +
466 454 # Run the script
467 455 self.run()
468 456  
... ... @@ -478,8 +466,5 @@ if __name__ == '__main__':
478 466 # Create instance of application
479 467 app = cstsdist(sys.argv)
480 468  
481   - # Open logfile
482   - app._logFileOpen()
483   -
484 469 # Execute application
485 470 app.execute()
... ...
cscripts/cstsmapmerge.par
  View file @31f23bc
... ... @@ -7,9 +7,9 @@
7 7 # File: cstsmapmerge.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
13 13 # #
14 14 # Purpose: Parameter file for the cstsmapmerge script #
15 15 # #
... ...
cscripts/cstsmapmerge.py
  View file @31f23bc
... ... @@ -33,7 +33,7 @@ class cstsmapmerge(ctools.cscript):
33 33 Merge Test Statistic maps.
34 34 """
35 35  
36   - # Constructors and destructors
  36 + # Constructor
37 37 def __init__(self, *argv):
38 38 """
39 39 Constructor.
... ... @@ -52,28 +52,11 @@ class cstsmapmerge(ctools.cscript):
52 52 self._merged_files = []
53 53 self._overwrite = True
54 54 self._delete = False
55   - self._outmap = "NONE"
56   - self._read_ahead = False
57   -
58   - # Initialise application
59   - if len(argv) == 0:
60   - ctools.cscript.__init__(self, self._name, self._version)
61   - elif len(argv) == 1:
62   - ctools.cscript.__init__(self, self._name, self._version, *argv)
63   - else:
64   - raise TypeError("Invalid number of arguments given.")
65   -
66   - # Set logger properties
67   - self._log_header()
68   - self._log.date(True)
69 55  
70   - # Return
71   - return
  56 + # Initialise application by calling the appropriate class
  57 + # constructor.
  58 + self._init_cscript(argv)
72 59  
73   - def __del__(self):
74   - """
75   - Destructor.
76   - """
77 60 # Return
78 61 return
79 62  
... ... @@ -120,8 +103,8 @@ class cstsmapmerge(ctools.cscript):
120 103 self._delete = self["delete"].boolean()
121 104  
122 105 # Read ahead output filename
123   - if self._read_ahead:
124   - self._outmap = self["outmap"].filename()
  106 + if self._read_ahead():
  107 + self["outmap"].filename()
125 108  
126 109 # Write input parameters into logger
127 110 if self._logTerse():
... ... @@ -391,12 +374,12 @@ class cstsmapmerge(ctools.cscript):
391 374 self._log.header1("Save TS map")
392 375  
393 376 # Get output filename in case it was not read ahead
394   - self._outmap = self["outmap"].filename()
  377 + outmap = self["outmap"].filename()
395 378  
396 379 # Log filename
397 380 if self._logTerse():
398 381 self._log(gammalib.parformat("TS map file"))
399   - self._log(self._outmap.url())
  382 + self._log(outmap.url())
400 383 self._log("\n")
401 384  
402 385 # Create FITS file
... ... @@ -426,7 +409,7 @@ class cstsmapmerge(ctools.cscript):
426 409 fits[fits.size()-1].extname("STATUS MAP")
427 410  
428 411 # Save FITS file
429   - fits.saveto(self._outmap, self._clobber())
  412 + fits.saveto(outmap, self._clobber())
430 413  
431 414 # Delete TS input maps if requested
432 415 if self._delete:
... ... @@ -444,8 +427,11 @@ class cstsmapmerge(ctools.cscript):
444 427 """
445 428 Execute the script.
446 429 """
447   - # Set read ahead flag
448   - self._read_ahead = True
  430 + # Open logfile
  431 + self._logFileOpen()
  432 +
  433 + # Read ahead output parameters
  434 + self._read_ahead(True)
449 435  
450 436 # Run the script
451 437 self.run()
... ... @@ -465,8 +451,5 @@ if __name__ == '__main__':
465 451 # Create instance of application
466 452 app = cstsmapmerge(sys.argv)
467 453  
468   - # Open logfile
469   - app._logFileOpen()
470   -
471 454 # Execute application
472 455 app.execute()
... ...
cscripts/obsutils.py
  View file @31f23bc
1 1 # ==========================================================================
2   -# This script provides a number of functions that are useful for handling
3   -# CTA observations.
  2 +# CTA observation handling support functions.
4 3 #
5   -# Copyright (C) 2011-2015 Juergen Knoedlseder
  4 +# Copyright (C) 2011-2016 Juergen Knoedlseder
6 5 #
7 6 # This program is free software: you can redistribute it and/or modify
8 7 # it under the terms of the GNU General Public License as published by
... ... @@ -31,21 +30,28 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
31 30 """
32 31 Simulate events for all observations in the container.
33 32  
34   - Parameters:
35   - obs - Observation container
36   - Keywords:
37   - log - Create log file(s)
38   - debug - Create console dump?
39   - edisp - Apply energy dispersion?
40   - seed - Seed value for simulations (default: 0)
41   - nbins - Number of energy bins (default: 0=unbinned)
42   - binsz - Pixel size for binned simulation (deg/pixel)
43   - npix - Number of pixels in X and Y for binned simulation
  33 + Args:
  34 + obs: Observation container without events
  35 +
  36 + Kwargs:
  37 + log: Create log file(s)
  38 + debug: Create console dump?
  39 + chatter: Chatter level
  40 + edisp: Apply energy dispersion?
  41 + seed: Seed value for simulations (default: 0)
  42 + nbins: Number of energy bins (default: 0=unbinned)
  43 + binsz: Pixel size for binned simulation (deg/pixel)
  44 + npix: Number of pixels in X and Y for binned simulation
  45 + proj: Projection for binned simulation
  46 + coord: Coordinate system for binned simulation
  47 +
  48 + Returns:
  49 + Observation container, filled with simulated events.
44 50 """
45 51  
46 52 # Allocate ctobssim application and set parameters
47 53 sim = ctools.ctobssim(obs)
48   - sim["seed"] = seed
  54 + sim["seed"] = seed
49 55 sim["edisp"] = edisp
50 56  
51 57 # Optionally open the log file
... ... @@ -85,16 +91,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
85 91  
86 92 # Allocate ctbin application and set parameters
87 93 bin = ctools.ctbin(sim.obs())
88   - bin["ebinalg"] = "LOG"
89   - bin["emin"] = emin
90   - bin["emax"] = emax
  94 + bin["ebinalg"] = "LOG"
  95 + bin["emin"] = emin
  96 + bin["emax"] = emax
91 97 bin["enumbins"] = nbins
92   - bin["usepnt"] = True # Use pointing for map centre
93   - bin["nxpix"] = npix
94   - bin["nypix"] = npix
95   - bin["binsz"] = binsz
  98 + bin["usepnt"] = True # Use pointing for map centre
  99 + bin["nxpix"] = npix
  100 + bin["nypix"] = npix
  101 + bin["binsz"] = binsz
96 102 bin["coordsys"] = coord
97   - bin["proj"] = proj
  103 + bin["proj"] = proj
98 104  
99 105 # Optionally open the log file
100 106 if log:
... ... @@ -350,28 +356,30 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
350 356 # ======================= #
351 357 def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
352 358 emin=0.1, emax=100.0, rad=5.0, \
353   - irf="South_50h", caldb="prod2", id="000000", instrument="CTA"):
  359 + irf="South_50h", caldb="prod2", id="000000"):
354 360 """
355   - Returns a single CTA observation containing an empty CTA event list.
356   - By looping over this function you can add CTA observations to the
357   - observation container.
358   -
359   - Parameters:
360   - pntdir - Pointing direction [GSkyDir]
361   - Keywords:
362   - tstart - Start time [seconds] (default: 0.0)
363   - duration - Duration of observation [seconds] (default: 1800.0)
364   - deadc - Deadtime correction factor (default: 0.95)
365   - emin - Minimum event energy [TeV] (default: 0.1)
366   - emax - Maximum event energy [TeV] (default: 100.0)
367   - rad - ROI radius used for analysis [deg] (default: 5.0)
368   - irf - Instrument response function (default: cta_dummy_irf)
369   - caldb - Calibration database path (default: "dummy")
370   - id - Run identifier (default: "000000")
371   - instrument - Intrument (default: "CTA")
  361 + Set a single CTA observation.
  362 +
  363 + The function sets a single CTA observation containing an empty CTA
  364 + event list. By looping over this function you can add CTA observations
  365 + to the observation container.
  366 +
  367 + Args:
  368 + pntdir: Pointing direction [GSkyDir]
  369 +
  370 + Kwargs:
  371 + tstart: Start time (seconds) (default: 0.0)
  372 + duration: Duration of observation (seconds) (default: 1800.0)
  373 + deadc: Deadtime correction factor (default: 0.95)
  374 + emin: Minimum event energy (TeV) (default: 0.1)
  375 + emax: Maximum event energy (TeV) (default: 100.0)
  376 + rad: ROI radius used for analysis (deg) (default: 5.0)
  377 + irf: Instrument response function (default: "South_50h")
  378 + caldb: Calibration database path (default: "prod2")
  379 + id: Run identifier (default: "000000")
372 380 """
373 381 # Allocate CTA observation
374   - obs_cta = gammalib.GCTAObservation(instrument)
  382 + obs_cta = gammalib.GCTAObservation()
375 383  
376 384 # Set calibration database
377 385 db = gammalib.GCaldb()
... ... @@ -397,7 +405,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
397 405 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
398 406  
399 407 # Set energy boundaries
400   - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
  408 + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"),
401 409 gammalib.GEnergy(emax, "TeV"))
402 410  
403 411 # Allocate event list
... ...
pyext/ctedispcube.i 0 → 100644
  View file @31f23bc
  1 +/***************************************************************************
  2 + * ctedispcube - Energy dispersion cube generation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Maria Haupt *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctedispcube.i
  23 + * @brief Energy dispersion cube generation tool definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "ctedispcube.hpp"
  30 +%}
  31 +
  32 +/***********************************************************************//**
  33 + * @class ctedispcube
  34 + *
  35 + * @brief Energy dispersion cube generation tool
  36 + ***************************************************************************/
  37 +class ctedispcube : public ctool {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + ctedispcube(void);
  42 + explicit ctedispcube(const GObservations& obs);
  43 + ctedispcube(int argc, char *argv[]);
  44 + ctedispcube(const ctedispcube& app);
  45 + virtual ~ctedispcube(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + void run(void);
  50 + void save(void);
  51 + const GCTACubeEdisp& edispcube(void) const;
  52 +};
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Energy dispersion cube generation tool Python extensions
  56 + ***************************************************************************/
  57 +%extend ctedispcube {
  58 + ctedispcube copy() {
  59 + return (*self);
  60 + }
  61 +}
... ...
pyext/ctool.i
  View file @31f23bc
... ... @@ -99,6 +99,7 @@ public:
99 99 %rename(_set_obs_bounds) set_obs_bounds(GObservations& obs);
100 100 %rename(_get_mean_pointing) get_mean_pointing(const GObservations& obs);
101 101 %rename(_get_current_rss) get_current_rss();
  102 + %rename(_get_obs_header) get_obs_header();
102 103  
103 104 // Protected methods
104 105 const bool& read_ahead(void) const;
... ... @@ -122,6 +123,7 @@ public:
122 123 void set_obs_bounds(GObservations& obs);
123 124 GSkyDir get_mean_pointing(const GObservations& obs);
124 125 size_t get_current_rss(void);
  126 + std::string get_obs_header(const GObservation* obs);
125 127 };
126 128  
127 129  
... ... @@ -152,7 +154,7 @@ public:
152 154  
153 155  
154 156 /***********************************************************************//**
155   - * @brief ctool base class Python extension
  157 + * @brief ctool base class Python extensions
156 158 ***************************************************************************/
157 159 %extend ctool {
158 160 void _read_ahead(const bool& flag) {
... ... @@ -160,3 +162,35 @@ public:
160 162 return;
161 163 }
162 164 }
  165 +
  166 +/***********************************************************************//**
  167 + * @brief cscript base class Python extensions
  168 + ***************************************************************************/
  169 +%pythoncode %{
  170 +
  171 +# Initialise observation container from constructor arguments. In case that
  172 +# an observation container is provided as argument, this container will be
  173 +# used to initialise the class.
  174 +def _set_input_obs(self, argv):
  175 + if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
  176 + obs = argv[0]
  177 + argv = argv[1:]
  178 + else:
  179 + obs = gammalib.GObservations()
  180 + return obs
  181 +cscript._set_input_obs = _set_input_obs
  182 +
  183 +# Initialise application by calling the appropriate class constructor.
  184 +def _init_cscript(self, argv):
  185 + if len(argv) == 0:
  186 + cscript.__init__(self, self._name, self._version)
  187 + elif len(argv) == 1:
  188 + cscript.__init__(self, self._name, self._version, *argv)
  189 + else:
  190 + raise TypeError("Invalid number of arguments given.")
  191 + # Set logger properties
  192 + self._log_header()
  193 + self._log.date(True)
  194 +cscript._init_cscript = _init_cscript
  195 +
  196 +%}
... ...
pyext/ctools/tools.i
  View file @31f23bc
1 1 /***************************************************************************
2 2 * ctools - SWIG file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -54,6 +54,7 @@
54 54 %include "ctskymap.i"
55 55 %include "ctexpcube.i"
56 56 %include "ctpsfcube.i"
  57 +%include "ctedispcube.i"
57 58 %include "ctbkgcube.i"
58 59 %include "ctcubemask.i"
59 60 %include "cttsmap.i"
... ... @@ -66,4 +67,4 @@
66 67 def test():
67 68 from ctools.tests import test_python_ctools
68 69 test_python_ctools.test(installed=True)
69   -%}
70 70 \ No newline at end of file
  71 +%}
... ...
src/ctbin/ctbin.cpp
  View file @31f23bc
... ... @@ -52,6 +52,8 @@ const GEnergy g_energy_margin(1.0e-12, &quot;TeV&quot;);
52 52  
53 53 /***********************************************************************//**
54 54 * @brief Void constructor
  55 + *
  56 + * Constructs an empty ctbin tool.
55 57 ***************************************************************************/
56 58 ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION)
57 59 {
... ... @@ -68,8 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION)
68 70 *
69 71 * param[in] obs Observation container.
70 72 *
71   - * This method creates an instance of the class by copying an existing
72   - * observations container.
  73 + * Constructs ctbin tool from an observation container.
73 74 ***************************************************************************/
74 75 ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION)
75 76 {
... ... @@ -90,6 +91,9 @@ ctbin::ctbin(const GObservations&amp; obs) : ctool(CTBIN_NAME, CTBIN_VERSION)
90 91 *
91 92 * @param[in] argc Number of arguments in command line.
92 93 * @param[in] argv Array of command line arguments.
  94 + *
  95 + * Constructs ctbin tool using command line arguments for user parameter
  96 + * setting.
93 97 ***************************************************************************/
94 98 ctbin::ctbin(int argc, char *argv[]) :
95 99 ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv)
... ... @@ -106,6 +110,8 @@ ctbin::ctbin(int argc, char *argv[]) :
106 110 * @brief Copy constructor
107 111 *
108 112 * @param[in] app Application.
  113 + *
  114 + * Constructs ctbin tool from another ctbin instance.
109 115 ***************************************************************************/
110 116 ctbin::ctbin(const ctbin& app) : ctool(app)
111 117 {
... ... @@ -122,6 +128,8 @@ ctbin::ctbin(const ctbin&amp; app) : ctool(app)
122 128  
123 129 /***********************************************************************//**
124 130 * @brief Destructor
  131 + *
  132 + * Destructs ctbin tool.
125 133 ***************************************************************************/
126 134 ctbin::~ctbin(void)
127 135 {
... ... @@ -142,7 +150,10 @@ ctbin::~ctbin(void)
142 150 /***********************************************************************//**
143 151 * @brief Assignment operator
144 152 *
145   - * @param[in] app Application.
  153 + * @param[in] app ctbin tool.
  154 + * @return ctbin tool.
  155 + *
  156 + * Assigns ctbin tool.
146 157 ***************************************************************************/
147 158 ctbin& ctbin::operator=(const ctbin& app)
148 159 {
... ... @@ -175,7 +186,9 @@ ctbin&amp; ctbin::operator=(const ctbin&amp; app)
175 186 ==========================================================================*/
176 187  
177 188 /***********************************************************************//**
178   - * @brief Clear instance
  189 + * @brief Clear ctbin tool
  190 + *
  191 + * Clears ctbin tool.
179 192 ***************************************************************************/
180 193 void ctbin::clear(void)
181 194 {
... ... @@ -195,12 +208,10 @@ void ctbin::clear(void)
195 208  
196 209  
197 210 /***********************************************************************//**
198   - * @brief Bin the event data
  211 + * @brief Run the ctbin tool
199 212 *
200   - * This method loops over all observations found in the observation conatiner
201   - * and bins all events from the event list(s) into counts map(s). Note that
202   - * each event list is binned in a separate counts map, hence no summing of
203   - * events is done.
  213 + * Gets the user parameters and loops over all CTA observations in the
  214 + * observation container to bin the events into a single counts cube.
204 215 ***************************************************************************/
205 216 void ctbin::run(void)
206 217 {
... ... @@ -212,33 +223,17 @@ void ctbin::run(void)
212 223 // Get task parameters
213 224 get_parameters();
214 225  
215   - // Write parameters into logger
216   - if (logTerse()) {
217   - log_parameters();
218   - log << std::endl;
219   - }
220   -
221 226 // Write observation(s) into logger
222 227 if (logTerse()) {
223 228 log << std::endl;
224   - if (m_obs.size() > 1) {
225   - log.header1("Observations");
226   - }
227   - else {
228   - log.header1("Observation");
229   - }
  229 + log.header1(gammalib::number("Observation", m_obs.size()));
230 230 log << m_obs << std::endl;
231 231 }
232 232  
233 233 // Write header
234 234 if (logTerse()) {
235 235 log << std::endl;
236   - if (m_obs.size() > 1) {
237   - log.header1("Bin observations");
238   - }
239   - else {
240   - log.header1("Bin observation");
241   - }
  236 + log.header1(gammalib::number("Bin observation", m_obs.size()));
242 237 }
243 238  
244 239 // Loop over all observations in the container
... ... @@ -246,14 +241,7 @@ void ctbin::run(void)
246 241  
247 242 // Write header for observation
248 243 if (logTerse()) {
249   - std::string header = m_obs[i]->instrument() + " observation";
250   - if (m_obs[i]->name().length() > 1) {
251   - header += " \"" + m_obs[i]->name() + "\"";
252   - }
253   - if (m_obs[i]->id().length() > 1) {
254   - header += " (id=" + m_obs[i]->id() +")";
255   - }
256   - log.header3(header);
  244 + log.header3(get_obs_header(m_obs[i]));
257 245 }
258 246  
259 247 // Get CTA observation
... ... @@ -293,12 +281,7 @@ void ctbin::run(void)
293 281 // Write observation(s) into logger
294 282 if (logTerse()) {
295 283 log << std::endl;
296   - if (m_obs.size() > 1) {
297   - log.header1("Binned observations");
298   - }
299   - else {
300   - log.header1("Binned observation");
301   - }
  284 + log.header1(gammalib::number("Binned observation", m_obs.size()));
302 285 log << m_obs << std::endl;
303 286 }
304 287  
... ... @@ -339,8 +322,8 @@ void ctbin::save(void)
339 322  
340 323 // Log filename
341 324 if (logTerse()) {
342   - log << "Save counts cube into file \""+m_outcube+"\".";
343   - log << std::endl;
  325 + log << gammalib::parformat("Counts cube file");
  326 + log << m_outcube.url() << std::endl;
344 327 }
345 328  
346 329 // Save cube
... ... @@ -376,7 +359,8 @@ void ctbin::publish(const std::string&amp; name)
376 359  
377 360 // Log filename
378 361 if (logTerse()) {
379   - log << "Publish \""+user_name+"\" counts cube." << std::endl;
  362 + log << gammalib::parformat("Counts cube name");
  363 + log << user_name << std::endl;
380 364 }
381 365  
382 366 // Publish counts cube
... ... @@ -496,6 +480,12 @@ void ctbin::get_parameters(void)
496 480 m_outcube = (*this)["outcube"].filename();
497 481 }
498 482  
  483 + // Write parameters into logger
  484 + if (logTerse()) {
  485 + log_parameters();
  486 + log << std::endl;
  487 + }
  488 +
499 489 // Return
500 490 return;
501 491 }
... ... @@ -521,8 +511,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
521 511 const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs->events());
522 512 if (events == NULL) {
523 513 std::string msg = "CTA Observation does not contain an event "
524   - "list. Event list information is needed to "
525   - "fill the counts map.";
  514 + "list. An event list is needed to fill the "
  515 + "counts cube.";
526 516 throw GException::invalid_value(G_FILL_CUBE, msg);
527 517 }
528 518  
... ... @@ -534,7 +524,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
534 524  
535 525 // Check for RoI sanity
536 526 if (!roi.is_valid()) {
537   - std::string msg = "No RoI information found in input observation "
  527 + std::string msg = "No valid RoI found in input observation "
538 528 "\""+obs->name()+"\". Run ctselect to specify "
539 529 "an RoI for this observation before running "
540 530 "ctbin.";
... ... @@ -624,10 +614,12 @@ void ctbin::fill_cube(GCTAObservation* obs)
624 614 * @brief Create output observation container.
625 615 *
626 616 * Creates an output observation container that combines all input CTA
627   - * observation into a single cube-style observation. All non CTA observations
628   - * present in the observation container are kept. The method furthermore
629   - * conserves any response information in case that a single CTA observation
630   - * is provided. This supports the original binned analysis.
  617 + * observation into a single stacked observation. All non-CTA observations
  618 + * and all binned CTA observations that were present in the observation
  619 + * container are append to the observation container so that they can
  620 + * be used by other tools. The method furthermore conserves any response
  621 + * information in case that a single CTA observation is provided to support
  622 + * binned analysis.
631 623 ***************************************************************************/
632 624 void ctbin::obs_cube(void)
633 625 {
... ... @@ -641,26 +633,37 @@ void ctbin::obs_cube(void)
641 633 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
642 634 if (obs != NULL) {
643 635  
644   - // Only change the event type if we had an unbinned observation
  636 + // Change the event type if we had an unbinned observation
645 637 if (obs->eventtype() == "EventList") {
646 638  
  639 + // Assign cube to the observation
647 640 obs->events(this->cube());
648   - obs->eventfile("");
649   - } // endif: observation was unbinned
650   - else { // input observation was binned and skipped
651 641  
652   - // Create new and empty cube
653   - GCTAEventCube cube = GCTAEventCube(m_cube, m_ebounds, obs->gti());
  642 + }
  643 +
  644 + // ... otherwise the input observation was binned and hence
  645 + // skipped. In that case we simply append an empty counts
  646 + // cube
  647 + else {
654 648  
655   - // Assign empty cube to have a new binned observation
  649 + // Create empty counts cube
  650 + GCTAEventCube cube(m_cube, m_ebounds, obs->gti());
  651 +
  652 + // Assign empty cube to observation
656 653 obs->events(cube);
  654 +
657 655 }
658 656  
  657 + // Reset file name
  658 + obs->eventfile("");
  659 +
659 660 } // endif: obervation was valid
660 661  
661 662 } // endif: we only had one observation in the container
662 663  
663   - // ... otherwise put a single CTA observation in container
  664 + // ... otherwise put a single CTA observation in container and
  665 + // append all observations that have not been used in the binning
  666 + // to the container
664 667 else {
665 668  
666 669 // Allocate observation container
... ... @@ -672,10 +675,25 @@ void ctbin::obs_cube(void)
672 675 // Attach event cube to CTA observation
673 676 obs.events(this->cube());
674 677  
675   - // Set map centre as pointing
676   - GSkyPixel pixel(0.5*double(m_cube.nx()), 0.5*double(m_cube.ny()));
677   - GSkyDir centre = m_cube.pix2dir(pixel);
678   - GCTAPointing pointing(centre);
  678 + // Compute average pointing direction for all CTA event lists
  679 + double ra = 0.0;
  680 + double dec = 0.0;
  681 + double number = 0.0;
  682 + for (int i = 0; i < m_obs.size(); ++i) {
  683 + GCTAObservation* cta = dynamic_cast<GCTAObservation*>(m_obs[i]);
  684 + if ((cta != NULL) && (cta->eventtype() == "EventList")) {
  685 + ra += cta->pointing().dir().ra();
  686 + dec += cta->pointing().dir().dec();
  687 + number += 1.0;
  688 + }
  689 + }
  690 + if (number > 0.0) {
  691 + ra /= number;
  692 + dec /= number;
  693 + }
  694 + GSkyDir dir;
  695 + dir.radec(ra, dec);
  696 + GCTAPointing pointing(dir);
679 697  
680 698 // Compute deadtime correction
681 699 double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0;
... ... @@ -683,8 +701,8 @@ void ctbin::obs_cube(void)
683 701 // Set CTA observation attributes
684 702 obs.pointing(pointing);
685 703 obs.obs_id(0);
686   - obs.ra_obj(centre.ra_deg()); //!< Dummy
687   - obs.dec_obj(centre.dec_deg()); //!< Dummy
  704 + obs.ra_obj(dir.ra_deg()); //!< Dummy
  705 + obs.dec_obj(dir.dec_deg()); //!< Dummy
688 706 obs.ontime(m_ontime);
689 707 obs.livetime(m_livetime);
690 708 obs.deadc(deadc);
... ... @@ -701,6 +719,9 @@ void ctbin::obs_cube(void)
701 719 if (obs == NULL) {
702 720 container.append(*m_obs[i]);
703 721 }
  722 + else if (obs->eventtype() != "EventList") {
  723 + container.append(*m_obs[i]);
  724 + }
704 725 }
705 726  
706 727 // Set observation container
... ...
src/ctbin/ctbin.par
  View file @31f23bc
... ... @@ -7,27 +7,27 @@
7 7 # File: ctbin.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Jurgen Knodlseder #
  10 +# Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
12 12 # IRAP #
13 13 # #
14   -# Purpose: Parameter file for the ctbin application #
  14 +# Purpose: Parameter file for the ctbin tool #
15 15 # #
16 16 #############################################################################
17 17 #
18 18 # General parameters
19 19 #===================
20   -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
21   -outcube, f, a, "cntcube.fits",,, "Output counts cube file"
  20 +inobs, f, a, events.fits,,, "Input event list or observation definition XML file"
  21 +outcube, f, a, cntcube.fits,,, "Output counts cube file"
22 22  
23 23 #
24 24 # Energy binning parameters
25 25 #==========================
26   -ebinalg, s, a, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins"
  26 +ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins"
27 27 emin, r, a, 0.1,,, "Start value for first energy bin in TeV"
28 28 emax, r, a, 100.0,,, "Stop value for last energy bin in TeV"
29   -enumbins, i, a, 20,,, "Number of energy bins"
30   -ebinfile, f, a, "NONE",,, "Name of the file containing the energy bin definition"
  29 +enumbins, i, a, 20,1,200, "Number of energy bins"
  30 +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
31 31  
32 32 #
33 33 # Spatial binning parameters
... ... @@ -36,10 +36,10 @@ usepnt, b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
36 36 nxpix, i, a, 200,,, "Size of the X axis in pixels"
37 37 nypix, i, a, 200,,, "Size of the Y axis in pixels"
38 38 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
39   -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  39 +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
40 40 xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
41 41 yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
42   -proj, s, a, "CAR",AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
  42 +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
43 43  
44 44 #
45 45 # Standard parameters
... ... @@ -48,5 +48,5 @@ publish, b, h, no,,, &quot;Publish counts cube on VO Hub?&quot;
48 48 chatter, i, h, 2,0,4, "Chattiness of output"
49 49 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
50 50 debug, b, h, no,,, "Debugging mode activated"
51   -mode, s, h, "ql",,, "Mode of automatic parameters"
52   -logfile, f, h, "ctbin.log",,, "Log filename"
  51 +mode, s, h, ql,,, "Mode of automatic parameters"
  52 +logfile, f, h, ctbin.log,,, "Log filename"
... ...
src/ctobssim/ctobssim.cpp
  View file @31f23bc
... ... @@ -60,6 +60,8 @@ const double g_roi_margin = 0.5; //!&lt; Simulation radius margin (degrees)
60 60  
61 61 /***********************************************************************//**
62 62 * @brief Void constructor
  63 + *
  64 + * Constructs an empty ctobssim tool.
63 65 ***************************************************************************/
64 66 ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
65 67 {
... ... @@ -76,7 +78,7 @@ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
76 78 *
77 79 * @param[in] obs Observation container.
78 80 *
79   - * Constructs application by initialising from an observation container.
  81 + * Constructs ctobssim tool from an observation container.
80 82 ***************************************************************************/
81 83 ctobssim::ctobssim(const GObservations& obs) :
82 84 ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
... ... @@ -98,7 +100,8 @@ ctobssim::ctobssim(const GObservations&amp; obs) :
98 100 * @param[in] argc Number of arguments in command line.
99 101 * @param[in] argv Array of command line arguments.
100 102 *
101   - * Constructs application using command line arguments for parameter setting.
  103 + * Constructs ctobssim tool using command line arguments for user parameter
  104 + * setting.
102 105 ***************************************************************************/
103 106 ctobssim::ctobssim(int argc, char *argv[]) :
104 107 ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
... ... @@ -115,6 +118,8 @@ ctobssim::ctobssim(int argc, char *argv[]) :
115 118 * @brief Copy constructor
116 119 *
117 120 * @param[in] app Application.
  121 + *
  122 + * Constructs ctobssim tool from another ctobssim instance.
118 123 ***************************************************************************/
119 124 ctobssim::ctobssim(const ctobssim& app) : ctool(app)
120 125 {
... ... @@ -131,6 +136,8 @@ ctobssim::ctobssim(const ctobssim&amp; app) : ctool(app)
131 136  
132 137 /***********************************************************************//**
133 138 * @brief Destructor
  139 + *
  140 + * Destructs ctobssim tool.
134 141 ***************************************************************************/
135 142 ctobssim::~ctobssim(void)
136 143 {
... ... @@ -151,8 +158,10 @@ ctobssim::~ctobssim(void)
151 158 /***********************************************************************//**
152 159 * @brief Assignment operator
153 160 *
154   - * @param[in] app Application.
155   - * @return Application.
  161 + * @param[in] app ctobssim tool.
  162 + * @return ctobssim tool.
  163 + *
  164 + * Assigns ctobssim tool.
156 165 ***************************************************************************/
157 166 ctobssim& ctobssim::operator=(const ctobssim& app)
158 167 {
... ... @@ -185,7 +194,9 @@ ctobssim&amp; ctobssim::operator=(const ctobssim&amp; app)
185 194 ==========================================================================*/
186 195  
187 196 /***********************************************************************//**
188   - * @brief Clear application
  197 + * @brief Clear ctobssim tool
  198 + *
  199 + * Clears ctobssim tool.
189 200 ***************************************************************************/
190 201 void ctobssim::clear(void)
191 202 {
... ... @@ -205,10 +216,10 @@ void ctobssim::clear(void)
205 216  
206 217  
207 218 /***********************************************************************//**
208   - * @brief Simulate event data
  219 + * @brief Run the ctobssim tool.
209 220 *
210   - * This method runs the simulation. Results are not saved by this method.
211   - * Invoke "save" to save the results.
  221 + * Gets the user parameters, loops over all CTA observations in the
  222 + * observation container, and simulate events for each observation.
212 223 ***************************************************************************/
213 224 void ctobssim::run(void)
214 225 {
... ... @@ -286,24 +297,14 @@ void ctobssim::run(void)
286 297 // Write observation(s) into logger
287 298 if (logTerse()) {
288 299 log << std::endl;
289   - if (m_obs.size() > 1) {
290   - log.header1("Observations");
291   - }
292   - else {
293   - log.header1("Observation");
294   - }
  300 + log.header1(gammalib::number("Observation", m_obs.size()));
295 301 log << m_obs << std::endl;
296 302 }
297 303  
298 304 // Write header
299 305 if (logTerse()) {
300 306 log << std::endl;
301   - if (m_obs.size() > 1) {
302   - log.header1("Simulate observations");
303   - }
304   - else {
305   - log.header1("Simulate observation");
306   - }
  307 + log.header1(gammalib::number("Simulate observation", m_obs.size()));
307 308 }
308 309  
309 310 // From here on the code can be parallelized if OpenMP support
... ... @@ -334,14 +335,7 @@ void ctobssim::run(void)
334 335  
335 336 // Write header for observation
336 337 if (logTerse()) {
337   - std::string header = m_obs[i]->instrument() + " observation";
338   - if (m_obs[i]->name().length() > 1) {
339   - header += " \"" + m_obs[i]->name() + "\"";
340   - }
341   - if (m_obs[i]->id().length() > 1) {
342   - header += " (id=" + m_obs[i]->id() +")";
343   - }
344   - wrklog.header3(header);
  338 + wrklog.header3(get_obs_header(m_obs[i]));
345 339 }
346 340  
347 341 // Get pointer on CTA observation
... ... @@ -1455,8 +1449,8 @@ void ctobssim::save_fits(void)
1455 1449  
1456 1450 // Log filename
1457 1451 if (logTerse()) {
1458   - log << "Save event list into file \""+m_outevents+"\".";
1459   - log << std::endl;
  1452 + log << gammalib::parformat("Event list file");
  1453 + log << m_outevents << std::endl;
1460 1454 }
1461 1455  
1462 1456 // Save observation into FITS file
... ... @@ -1483,7 +1477,7 @@ void ctobssim::save_xml(void)
1483 1477 {
1484 1478 // Get output filename and prefix
1485 1479 m_outevents = (*this)["outevents"].filename();
1486   - m_prefix = (*this)["prefix"].string();
  1480 + m_prefix = (*this)["prefix"].string();
1487 1481  
1488 1482 // Issue warning if output filename has no .xml suffix
1489 1483 std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4));
... ... @@ -1513,15 +1507,16 @@ void ctobssim::save_xml(void)
1513 1507 if (obs->events()->size() != 0) {
1514 1508  
1515 1509 // Set event output file name
1516   - std::string outfile = m_prefix + gammalib::str(i) + ".fits";
  1510 + std::string outfile = m_prefix + gammalib::str(i) +
  1511 + ".fits";
1517 1512  
1518 1513 // Store output file name in observation
1519 1514 obs->eventfile(outfile);
1520 1515  
1521 1516 // Log filename
1522 1517 if (logTerse()) {
1523   - log << "Save event list into file \""+outfile+"\".";
1524   - log << std::endl;
  1518 + log << gammalib::parformat("Event list file");
  1519 + log << outfile << std::endl;
1525 1520 }
1526 1521  
1527 1522 // Save observation into FITS file
... ...
src/ctobssim/ctobssim.par
  View file @31f23bc
... ... @@ -7,7 +7,7 @@
7 7 # File: ctobssim.par #
8 8 # Version: 1.1.0 #
9 9 # #
10   -# Author: Jurgen Knodlseder #
  10 +# Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
12 12 # IRAP #
13 13 # #
... ... @@ -17,13 +17,13 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -inobs, f, h, "NONE",,, "Input event list or observation definition XML file"
21   -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
22   -caldb, s, a, "prod2",,, "Calibration database"
23   -irf, s, a, "South_0.5h",,, "Instrument response function"
  20 +inobs, f, h, NONE,,, "Input event list or observation definition XML file"
  21 +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file"
  22 +caldb, s, a, prod2,,, "Calibration database"
  23 +irf, s, a, South_0.5h,,, "Instrument response function"
24 24 edisp, b, h, no,,, "Apply energy dispersion?"
25   -outevents, f, a, "events.fits",,, "Output event data file or observation definition XML file"
26   -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition XML file"
  25 +outevents, f, a, events.fits,,, "Output event data file or observation definition XML file"
  26 +prefix, s, h, sim_events_,,, "Prefix for event lists in observation definition XML file"
27 27  
28 28 #
29 29 # Simulation parameters
... ... @@ -38,7 +38,7 @@ emin, r, a, 0.1,,, &quot;Lower energy limit (TeV)&quot;
38 38 emax, r, a, 100.0,,, "Upper energy limit (TeV)"
39 39 deadc, r, h, 0.95,0,1, "Average deadtime correction factor"
40 40 maxrate, r, h, 1.0e6,,, "Maximum photon rate"
41   -eslices, i, h, 10,1,100,"Number of energy slices"
  41 +eslices, i, h, 10,1,100, "Number of energy slices"
42 42  
43 43 #
44 44 # Standard parameters
... ... @@ -46,5 +46,5 @@ eslices, i, h, 10,1,100,&quot;Number of energy slices&quot;
46 46 chatter, i, h, 2,0,4, "Chattiness of output"
47 47 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
48 48 debug, b, h, no,,, "Debugging mode activated"
49   -mode, s, h, "ql",,, "Mode of automatic parameters"
50   -logfile, f, h, "ctobssim.log",,, "Log filename"
  49 +mode, s, h, ql,,, "Mode of automatic parameters"
  50 +logfile, f, h, ctobssim.log,,, "Log filename"
... ...
src/support/ctool.cpp
  View file @31f23bc
... ... @@ -1108,6 +1108,39 @@ size_t ctool::get_current_rss(void)
1108 1108  
1109 1109  
1110 1110 /***********************************************************************//**
  1111 + * @brief Return observation header string
  1112 + *
  1113 + * @param[in] obs Pointer to observation.
  1114 + * @return String containing observation information.
  1115 + *
  1116 + * Returns a string that contains observation information, including the
  1117 + * instrument name, the observation name, and the observation ID. The format
  1118 + * of the string is
  1119 + *
  1120 + * XXX observation "name" (id=YYYYY)
  1121 + *
  1122 + ***************************************************************************/
  1123 +std::string ctool::get_obs_header(const GObservation* obs)
  1124 +{
  1125 + // Initialise header string
  1126 + std::string header = obs->instrument() + " observation";
  1127 +
  1128 + // If observation name is not empty then add name
  1129 + if (!obs->name().empty()) {
  1130 + header += " \"" + obs->name() + "\"";
  1131 + }
  1132 +
  1133 + // If observation ID is not empty then add ID
  1134 + if (!obs->id().empty()) {
  1135 + header += " (id=" + obs->id() +")";
  1136 + }
  1137 +
  1138 + // Return header
  1139 + return header;
  1140 +}
  1141 +
  1142 +
  1143 +/***********************************************************************//**
1111 1144 * @brief Dumps help text in the console
1112 1145 ***************************************************************************/
1113 1146 void ctool::provide_help(void) const
... ...
src/support/ctool.hpp
  View file @31f23bc
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -90,17 +90,18 @@ protected:
90 90 void set_obs_bounds(GObservations& obs);
91 91 GSkyDir get_mean_pointing(const GObservations& obs);
92 92 size_t get_current_rss(void);
  93 + std::string get_obs_header(const GObservation* obs);
93 94  
94 95 // Protected members
95   - bool m_read_ahead; //!< Read ahead parameters
  96 + bool m_read_ahead; //!< Read ahead output parameters
96 97  
97 98 protected:
98 99 // Protected methods
99   - void provide_help(void) const;
  100 + void provide_help(void) const;
100 101  
101 102 // Protected members
102   - bool m_use_xml; //!< Use XML file instead of FITS file for observations
103   - GTimeReference m_cta_ref; //!< CTA time reference
  103 + bool m_use_xml; //!< Use XML file instead of FITS file for observations
  104 + GTimeReference m_cta_ref; //!< CTA time reference
104 105 };
105 106  
106 107  
... ...
test/test_ctbin.py
  View file @31f23bc
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This scripts performs unit tests for the ctbin tool.
  3 +# This scripts performs unit tests for the ctbin module.
4 4 #
5 5 # Copyright (C) 2014-2016 Juergen Knoedlseder
6 6 #
... ... @@ -22,13 +22,14 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ========================= #
26   -# Test class for ctbin tool #
27   -# ========================= #
  25 +# =========================== #
  26 +# Test class for ctbin module #
  27 +# =========================== #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30   - Test class for ctbin tool.
  30 + Test class for ctbin module.
31 31 """
  32 +
32 33 # Constructor
33 34 def __init__(self):
34 35 """
... ... @@ -38,7 +39,7 @@ class Test(gammalib.GPythonTestSuite):
38 39 gammalib.GPythonTestSuite.__init__(self)
39 40  
40 41 # Set members
41   - self.events_name = "data/crab_events.fits"
  42 + self._events_name = "data/crab_events.fits"
42 43  
43 44 # Return
44 45 return
... ... @@ -52,19 +53,19 @@ class Test(gammalib.GPythonTestSuite):
52 53 self.name("ctbin")
53 54  
54 55 # Append tests
55   - self.append(self.test_functional, "Test ctbin functionality")
  56 + self.append(self._test_ctbin, "Test ctbin")
56 57  
57 58 # Return
58 59 return
59 60  
60   - # Test ctbin functionnality
61   - def test_functional(self):
  61 + # Test ctbin
  62 + def _test_ctbin(self):
62 63 """
63   - Test ctbin functionnality.
  64 + Test ctbin.
64 65 """
65 66 # Set-up ctbin
66 67 bin = ctools.ctbin()
67   - bin["inobs"] = self.events_name
  68 + bin["inobs"] = self._events_name
68 69 bin["outcube"] = "cntmap.fits"
69 70 bin["ebinalg"] = "LOG"
70 71 bin["emin"] = 0.1
... ... @@ -78,27 +79,154 @@ class Test(gammalib.GPythonTestSuite):
78 79 bin["xref"] = 83.63
79 80 bin["yref"] = 22.01
80 81  
81   - # Run tool
82   - self.test_try("Run ctbin")
83   - try:
84   - bin.run()
85   - self.test_try_success()
86   - #except Exception as e:
87   - # msg = "Exception occured in ctbin: %s." % (e,)
88   - except:
89   - msg = "Exception occured in ctbin."
90   - self.test_try_failure(msg)
  82 + # Run ctbin tool
  83 + bin.run()
  84 +
  85 + # Check content of observation and cube
  86 + self._test_observation(bin, 5542)
  87 + self._test_cube(bin.cube(), 5542)
  88 +
  89 + # Test copy constructor
  90 + cpy_bin = bin.copy()
  91 +
  92 + # Check content of observation and cube
  93 + self._test_observation(cpy_bin, 5542)
  94 + self._test_cube(cpy_bin.cube(), 5542)
  95 +
  96 + # Run copy of ctbin tool again
  97 + cpy_bin.run()
  98 +
  99 + # Check content of observation and cube. We expect now an empty
  100 + # event cube as on input the observation is binned, and any binned
  101 + # observation will be skipped, hence the counts cube should be
  102 + # empty.
  103 + self._test_observation(cpy_bin, 0)
  104 + self._test_cube(cpy_bin.cube(), 0)
  105 +
  106 + # Clear and run copy of ctbin again
  107 + #cpy_bin.clear()
  108 + #cpy_bin.run()
  109 +#RuntimeError: *** ERROR in GXmlAttribute::value(std::string): Invalid XML attribute value: *** ERROR in GPythonTestSuite::test: <type 'exceptions.RuntimeError'>
  110 +#*** ERROR in GApplicationPars::operator[](std::string&): Invalid argument. Parameter "debug" has not been found in parameter file.
  111 +#Please specify a valid parameter name.
91 112  
92 113 # Save counts cube
93   - self.test_try("Save counts cube")
  114 + bin.save()
  115 +
  116 + # Load counts cube and check content.
  117 + evt = gammalib.GCTAEventCube("cntmap.fits")
  118 + self._test_cube(evt, 5542)
  119 +
  120 + # Prepare observation container for stacked analysis
  121 + cta = gammalib.GCTAObservation(self._events_name)
  122 + obs = gammalib.GObservations()
  123 + cta.id("0001")
  124 + obs.append(cta)
  125 + cta.id("0002")
  126 + obs.append(cta)
  127 + cta.id("0003")
  128 + obs.append(cta)
  129 +
  130 + # Set-up ctbin using an observation container
  131 + bin = ctools.ctbin(obs)
  132 + bin["outcube"] = "cntmap.fits"
  133 + bin["ebinalg"] = "LOG"
  134 + bin["emin"] = 0.1
  135 + bin["emax"] = 100.0
  136 + bin["enumbins"] = 20
  137 + bin["nxpix"] = 200
  138 + bin["nypix"] = 200
  139 + bin["binsz"] = 0.02
  140 + bin["coordsys"] = "CEL"
  141 + bin["proj"] = "CAR"
  142 + bin["xref"] = 83.63
  143 + bin["yref"] = 22.01
  144 +
  145 + # Run ctbin tool
  146 + bin.run()
  147 +
  148 + # Check content of observation and cube (need multiplier=3 since
  149 + # three identical observations have been appended)
  150 + self._test_observation(bin, 5542, multiplier=3)
  151 + self._test_cube(bin.cube(), 5542, multiplier=3)
  152 +
  153 + # Set-up ctbin with invalid event file
  154 + bin = ctools.ctbin()
  155 + bin["inobs"] = "event_file_that_does_not_exist.fits"
  156 + bin["outcube"] = "cntmap.fits"
  157 + bin["ebinalg"] = "LOG"
  158 + bin["emin"] = 0.1
  159 + bin["emax"] = 100.0
  160 + bin["enumbins"] = 20
  161 + bin["nxpix"] = 200
  162 + bin["nypix"] = 200
  163 + bin["binsz"] = 0.02
  164 + bin["coordsys"] = "CEL"
  165 + bin["proj"] = "CAR"
  166 + bin["xref"] = 83.63
  167 + bin["yref"] = 22.01
  168 +
  169 + # Run ctbin tool
  170 + self.test_try("Run ctbin with invalid event file")
94 171 try:
95   - bin.save()
96   - self.test_try_success()
97   - #except Exception as e:
98   - # msg = "Exception occured in saving counts cube: %s." % (e,)
  172 + bin.run()
  173 + self.test_try_failure()
99 174 except:
100   - msg = "Exception occured in saving counts cube."
101   - self.test_try_failure(msg)
  175 + self.test_try_success()
  176 +
  177 + # Return
  178 + return
  179 +
  180 + # Check observation
  181 + def _test_observation(self, ctbin, nevents, multiplier=1):
  182 + """
  183 + Test content of an observation.
  184 +
  185 + Args:
  186 + ctbin: ctbin instance
  187 + nevents: Expected number of events
  188 +
  189 + Kwargs:
  190 + multiplier: Observation multiplier
  191 + """
  192 + # Test observation container
  193 + obs = gammalib.GCTAObservation(ctbin.obs()[0])
  194 + pnt = obs.pointing()
  195 + self.test_value(ctbin.obs().size(), 1,
  196 + "There is one observation")
  197 + self.test_assert(obs.instrument() == "CTA",
  198 + "Observation is CTA observation")
  199 + self.test_value(obs.ontime(), 1800.0*multiplier, 1.0e-6,
  200 + "Ontime is 1800 sec")
  201 + self.test_value(obs.livetime(), 1710.0*multiplier, 1.0e-6,
  202 + "Livetime is 1710 sec")
  203 + self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6,
  204 + "Pointing Right Ascension is 83.63 deg")
  205 + self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6,
  206 + "Pointing Declination is 22.01 deg")
  207 +
  208 + # Test event cube
  209 + self._test_cube(obs.events(), nevents, multiplier=multiplier)
  210 +
  211 + # Return
  212 + return
  213 +
  214 + # Check event cube
  215 + def _test_cube(self, cube, nevents, multiplier=1):
  216 + """
  217 + Test content of event cube."
  218 +
  219 + Args:
  220 + cube: Event cube
  221 + nevents: Expected number of events.
  222 +
  223 + Kwargs:
  224 + multiplier: Event number multiplier
  225 + """
  226 + # Test event cube
  227 + self.test_value(cube.size(), 800000, "800000 event bins")
  228 + self.test_value(cube.number(), nevents*multiplier,
  229 + str(nevents*multiplier)+" events")
102 230  
103 231 # Return
104 232 return
... ...
test/test_ctobssim.py
  View file @31f23bc
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This scripts performs unit tests for the ctobssim tool.
  3 +# This scripts performs unit tests for the ctobssim module.
4 4 #
5 5 # Copyright (C) 2014-2016 Juergen Knoedlseder
6 6 #
... ... @@ -20,15 +20,17 @@
20 20 # ==========================================================================
21 21 import gammalib
22 22 import ctools
  23 +from cscripts import obsutils
23 24  
24 25  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  26 +# ============================== #
  27 +# Test class for ctobssim module #
  28 +# ============================== #
28 29 class Test(gammalib.GPythonTestSuite):
29 30 """
30   - Test class for ctobssim tool.
  31 + Test class for ctobssim module.
31 32 """
  33 +
32 34 # Constructor
33 35 def __init__(self):
34 36 """
... ... @@ -38,9 +40,9 @@ class Test(gammalib.GPythonTestSuite):
38 40 gammalib.GPythonTestSuite.__init__(self)
39 41  
40 42 # Set members
41   - self.model_name = "data/crab.xml"
42   - self.caldb = "irf"
43   - self.irf = "cta_dummy_irf"
  43 + self._model_name = "data/crab.xml"
  44 + self._caldb = "irf"
  45 + self._irf = "cta_dummy_irf"
44 46  
45 47 # Return
46 48 return
... ... @@ -54,23 +56,23 @@ class Test(gammalib.GPythonTestSuite):
54 56 self.name("ctobssim")
55 57  
56 58 # Append tests
57   - self.append(self.test_functional, "Test ctobssim functionality")
58   - self.append(self.test_container, "Test ctobssim on observation container")
  59 + self.append(self._test_ctobssim, "Test ctobssim")
  60 + #self.append(self._test_container, "Test ctobssim on observation container")
59 61  
60 62 # Return
61 63 return
62 64  
63   - # Test ctobssim functionnality
64   - def test_functional(self):
  65 + # Test ctobssim
  66 + def _test_ctobssim(self):
65 67 """
66   - Test ctobssim functionnality.
  68 + Test ctobssim.
67 69 """
68 70 # Set-up ctobssim
69 71 sim = ctools.ctobssim()
70   - sim["inmodel"] = self.model_name
  72 + sim["inmodel"] = self._model_name
71 73 sim["outevents"] = "events.fits"
72   - sim["caldb"] = self.caldb
73   - sim["irf"] = self.irf
  74 + sim["caldb"] = self._caldb
  75 + sim["irf"] = self._irf
74 76 sim["ra"] = 83.63
75 77 sim["dec"] = 22.01
76 78 sim["rad"] = 10.0
... ... @@ -80,154 +82,123 @@ class Test(gammalib.GPythonTestSuite):
80 82 sim["emax"] = 100.0
81 83  
82 84 # Run tool
83   - self.test_try("Run ctobssim")
84   - try:
85   - sim.run()
86   - self.test_try_success()
87   - #except Exception as e:
88   - # msg = "Exception occured in ctobssim: %s." % (e,)
89   - except:
90   - msg = "Exception occured in ctobssim."
91   - self.test_try_failure(msg)
  85 + sim.run()
92 86  
93   - # Retrieve observation and check content
94   - obs = gammalib.GCTAObservation(sim.obs()[0])
95   - pnt = obs.pointing()
96   - evt = obs.events()
97   - self.test_assert(obs.instrument() == "CTA", "Observation not a CTA observation")
98   - self.test_value(sim.obs().size(), 1, "There is not a single observation")
99   - self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is not 1800 sec")
100   - self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is not 1710 sec")
101   - self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg")
102   - self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg")
103   - self.test_value(evt.size(), 4119, "Number of events is not 4119")
  87 + # Check content of observation
  88 + self._test_observation(sim)
  89 + self._test_list(sim.obs()[0].events(), 4119)
104 90  
105 91 # Save events
106   - self.test_try("Save events")
107   - try:
108   - sim.save()
109   - self.test_try_success()
110   - #except Exception as e:
111   - # msg = "Exception occured in saving events: %s." % (e,)
112   - except:
113   - msg = "Exception occured in saving events."
114   - self.test_try_failure(msg)
  92 + sim.save()
  93 +
  94 + # Load counts cube and check content.
  95 + evt = gammalib.GCTAEventList("events.fits")
  96 + self._test_list(evt, 4119)
115 97  
116   - # Test ctobssim on observation container
117   - def test_container(self):
118   - """
119   - Test ctobssim on observation container.
120   - """
121 98 # Set-up observation container
122   - obs = self.set_obs(4)
  99 + pnts = [{'ra': 83.63, 'dec': 21.01},
  100 + {'ra': 84.63, 'dec': 22.01},
  101 + {'ra': 83.63, 'dec': 23.01},
  102 + {'ra': 82.63, 'dec': 22.01}]
  103 + obs = obsutils.set_obs_list(pnts, caldb=self._caldb, irf=self._irf)
123 104  
124   - # Set-up ctobssim
  105 + # Set-up ctobssim from observation container
125 106 sim = ctools.ctobssim(obs)
126 107 sim["outevents"] = "sim_events.xml"
127   - sim["inmodel"] = self.model_name
  108 + sim["inmodel"] = self._model_name
  109 +
128 110 # Run tool
129   - self.test_try("Run ctobssim")
130   - try:
131   - sim.run()
132   - self.test_try_success()
133   - #except Exception as e:
134   - # msg = "Exception occured in ctobssim: %s." % (e,)
135   - except:
136   - msg = "Exception occured in ctobssim."
137   - self.test_try_failure(msg)
  111 + sim.run()
138 112  
139 113 # Retrieve observation and check content
140   - self.test_value(sim.obs().size(), 4, "There are not 4 observations")
  114 + self._test_observation(sim, nobs=4, pnts=pnts)
  115 + self._test_list(sim.obs()[0].events(), 4063)
  116 + self._test_list(sim.obs()[1].events(), 4059)
  117 + self._test_list(sim.obs()[2].events(), 4010)
  118 + self._test_list(sim.obs()[3].events(), 4138)
141 119  
142 120 # Save events
143   - self.test_try("Save events")
  121 + sim.save()
  122 +
  123 + # Load events
  124 + obs = gammalib.GObservations("sim_events.xml")
  125 +
  126 + # Retrieve observation and check content
  127 + self._test_list(obs[0].events(), 4063)
  128 + self._test_list(obs[1].events(), 4059)
  129 + self._test_list(obs[2].events(), 4010)
  130 + self._test_list(obs[3].events(), 4138)
  131 +
  132 + # Set-up ctobssim with invalid event file
  133 + sim = ctools.ctobssim()
  134 + sim["inmodel"] = "model_file_that_does_not_exist.xml"
  135 + sim["outevents"] = "events.fits"
  136 + sim["caldb"] = self._caldb
  137 + sim["irf"] = self._irf
  138 + sim["ra"] = 83.63
  139 + sim["dec"] = 22.01
  140 + sim["rad"] = 10.0
  141 + sim["tmin"] = 0.0
  142 + sim["tmax"] = 1800.0
  143 + sim["emin"] = 0.1
  144 + sim["emax"] = 100.0
  145 +
  146 + # Run ctbin tool
  147 + self.test_try("Run ctobssim with invalid model file")
144 148 try:
145   - sim.save()
146   - self.test_try_success()
147   - #except Exception as e:
148   - # msg = "Exception occured in saving events: %s." % (e,)
  149 + sim.run()
  150 + self.test_try_failure()
149 151 except:
150   - msg = "Exception occured in saving events."
151   - self.test_try_failure(msg)
  152 + self.test_try_success()
  153 +
  154 + # Return
  155 + return
152 156  
153   - # Setup observation container
154   - def set_obs(self, number):
  157 + # Check observation
  158 + def _test_observation(self, ctobssim, nobs=1,
  159 + pnts=[{'ra': 83.63, 'dec': 22.01}]):
155 160 """
156   - Setup observation container.
  161 + Test content of an observation.
  162 +
  163 + Args:
  164 + ctobssim: ctobssim instance
  165 +
  166 + Kwargs:
  167 + nobs: Number of observations
  168 + pnts: List of pointing dictionaries
157 169 """
158   - # Initialise empty observation container
159   - obs = gammalib.GObservations()
160   -
161   - # Initialise first time and identifier
162   - tstart = 0.0
163   - offset = -float(number)/2.0
164   -
165   - # Loop over number of observations
166   - for i in range(number):
167   - id = "%6.6d" % i
168   - obs_cta = self.set_one_obs(id, offset)
169   - offset += 1.0
170   - tstart += 1830.0
171   - obs.append(obs_cta)
172   -
173   - # Return observation container
174   - return obs
175   -
176   - # Setup one observation
177   - def set_one_obs(self, id, offset, \
178   - tstart=0.0, duration=1800.0, deadc=0.95, \
179   - emin=0.1, emax=100.0, rad=5.0, \
180   - irf="cta_dummy_irf", caldb="dummy"):
  170 + # Test observation container
  171 + self.test_value(ctobssim.obs().size(), nobs,
  172 + "There is one observation")
  173 + for i in range(ctobssim.obs().size()):
  174 + obs = gammalib.GCTAObservation(ctobssim.obs()[i])
  175 + pnt = obs.pointing()
  176 + self.test_assert(obs.instrument() == "CTA",
  177 + "Observation is CTA observation")
  178 + self.test_value(obs.ontime(), 1800.0, 1.0e-6,
  179 + "Ontime is 1800 sec")
  180 + self.test_value(obs.livetime(), 1710.0, 1.0e-6,
  181 + "Livetime is 1710 sec")
  182 + self.test_value(pnt.dir().ra_deg(), pnts[i]['ra'], 1.0e-6,
  183 + "Pointing Right Ascension is "+str(pnts[i]['ra'])+" deg")
  184 + self.test_value(pnt.dir().dec_deg(), pnts[i]['dec'], 1.0e-6,
  185 + "Pointing Declination is "+str(pnts[i]['dec'])+" deg")
  186 +
  187 + # Return
  188 + return
  189 +
  190 + # Check event list
  191 + def _test_list(self, list, nevents):
181 192 """
182   - Setup one observation for test purposes.
  193 + Test content of event list."
  194 +
  195 + Args:
  196 + list: Event list
  197 + nevents: Expected number of events.
183 198 """
184   - # Allocate CTA observation
185   - obs_cta = gammalib.GCTAObservation()
186   -
187   - # Set calibration database
188   - db = gammalib.GCaldb()
189   - if (gammalib.dir_exists(self.caldb)):
190   - db.rootdir(self.caldb)
191   - else:
192   - db.open("cta", self.caldb)
193   -
194   - # Set pointing direction
195   - pntdir = gammalib.GSkyDir()
196   - pntdir.radec_deg(83.63, 22.01+offset)
197   - pnt = gammalib.GCTAPointing()
198   - pnt.dir(pntdir)
199   - obs_cta.pointing(pnt)
200   -
201   - # Set ROI
202   - roi = gammalib.GCTARoi()
203   - instdir = gammalib.GCTAInstDir()
204   - instdir.dir(pntdir)
205   - roi.centre(instdir)
206   - roi.radius(rad)
207   -
208   - # Set GTI
209   - gti = gammalib.GGti()
210   - gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
211   -
212   - # Set energy boundaries
213   - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
214   - gammalib.GEnergy(emax, "TeV"))
215   -
216   - # Allocate event list
217   - events = gammalib.GCTAEventList()
218   - events.roi(roi)
219   - events.gti(gti)
220   - events.ebounds(ebounds)
221   - obs_cta.events(events)
222   -
223   - # Set instrument response
224   - obs_cta.response(self.irf, db)
225   -
226   - # Set ontime, livetime, and deadtime correction factor
227   - obs_cta.ontime(duration)
228   - obs_cta.livetime(duration*deadc)
229   - obs_cta.deadc(deadc)
230   - obs_cta.id(id)
231   -
232   - # Return CTA observation
233   - return obs_cta
  199 + # Test event list
  200 + self.test_value(list.size(), nevents, str(nevents)+" elements")
  201 + self.test_value(list.number(), nevents, str(nevents)+" events")
  202 +
  203 + # Return
  204 + return
... ...
test/test_python_cscripts.py
  View file @31f23bc
... ... @@ -34,7 +34,7 @@ def test(installed=False):
34 34 Perform unit testing for Python interface.
35 35 """
36 36 # Allocate test suite container
37   - suites = gammalib.GTestSuites("cscripts unit testing")
  37 + suites = gammalib.GTestSuites("cscripts Python module unit testing")
38 38  
39 39 # Allocate test suites and append them to the container
40 40 suite_cscaldb = test_cscaldb.Test()
... ...
test/test_python_ctools.py
  View file @31f23bc
... ... @@ -48,7 +48,7 @@ def test(installed=False):
48 48 Perform unit testing for Python interface.
49 49 """
50 50 # Allocate test suite container
51   - suites = gammalib.GTestSuites("ctools unit testing")
  51 + suites = gammalib.GTestSuites("ctools Python module unit testing")
52 52  
53 53 # Allocate test suites and append them to the container
54 54 suite_ctobssim = test_ctobssim.Test()
... ...