Marianne Lemoine / ctools

Commit ea3e9e13ba4aa10722110554d1258e66608d691e

Authored by Michael Mayer
1 parent 82ecc3ad

Added csiactobs and csiactdata and created a unit test - removed cshessobs

Inline Side-by-side
Showing 9 changed files with 860 additions and 635 deletions
cscripts/Makefile.am
  View file @ea3e9e1
@@ -3,14 +3,14 @@ @@ -3,14 +3,14 @@
3 # Scripts to be distributed in $(pyexecdir)/cscripts 3 # Scripts to be distributed in $(pyexecdir)/cscripts
4 pythonscripts = __init__.py \ 4 pythonscripts = __init__.py \
5 cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \ 5 cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
6 - cshessobs.py csobsdef.py csspec.py cslightcrv.py csobsinfo.py \ 6 + csiactobs.py csiactdata.py csobsdef.py csspec.py cslightcrv.py csobsinfo.py \
7 csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \ 7 csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \
8 csobs2caldb.py \ 8 csobs2caldb.py \
9 obsutils.py 9 obsutils.py
10 10
11 # Scripts for which symbolic links will be created in bindir 11 # Scripts for which symbolic links will be created in bindir
12 linkscripts = cspull cssens cstsdist cscaldb csresmap \ 12 linkscripts = cspull cssens cstsdist cscaldb csresmap \
13 - cshessobs csobsdef csspec cslightcrv csobsinfo \ 13 + csiactobs csiactdata csobsdef csspec cslightcrv csobsinfo \
14 csinfo csmodelinfo csmodelmerge cstsmapmerge \ 14 csinfo csmodelinfo csmodelmerge cstsmapmerge \
15 csobs2caldb 15 csobs2caldb
16 16
cscripts/cshessobs.py deleted
View file @82ecc3a
1 -#!/usr/bin/env python  
2 -# ==========================================================================  
3 -# Generation of an H.E.S.S. observation definition file.  
4 -#  
5 -# Copyright (C) 2015 Michael Mayer  
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 -import gammalib  
22 -import ctools  
23 -import sys  
24 -import os  
25 -import glob  
26 -import math  
27 -  
28 -  
29 -# =============== #  
30 -# cshessobs class #  
31 -# =============== #  
32 -class cshessobs(ctools.cscript):  
33 - """  
34 - This class implements the creation of a observation xml file for HESS  
35 - data analysis. This class is dedicated for use inside the HESS Collaboration,  
36 - i.e. it can only be used if you have access to HESS data in FITS format.  
37 - The HESSFITS environment has to be set correctly to use this script  
38 - Please see the internal documentation for this purpose:  
39 - https://hess-confluence.desy.de/confluence/display/HESS/HESS+Open+Source+Tools+-+HESS+data+and+IRFs+in+FITS+format  
40 - """  
41 - def __init__(self, *argv):  
42 - """  
43 - Constructor.  
44 - """  
45 - # Set name  
46 - self.name = "cshessobs"  
47 - self.version = "0.2.0"  
48 - self.datapath = ""  
49 - self.m_erange = gammalib.GEbounds()  
50 -  
51 - # Check on existence of HESSFITS environment variable  
52 - try:  
53 - self.datapath = os.environ["HESSFITS"]  
54 - except KeyError:  
55 - raise RuntimeError("Environment variable \"HESSFITS\" not set. " \  
56 - "Please set the environment variable \"HESSFITS\" to use this script")  
57 -  
58 - # Initialise some members  
59 - self.inmodels = None  
60 - self.outobs = "obs.xml"  
61 -  
62 - # Make sure that parfile exists  
63 - file = self.parfile()  
64 -  
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)  
76 -  
77 - # Return  
78 - return  
79 -  
80 - def __del__(self):  
81 - """  
82 - Destructor.  
83 - """  
84 - # Return  
85 - return  
86 -  
87 - def parfile(self):  
88 - """  
89 - Check if parfile exists. If parfile does not exist then create a  
90 - default parfile. This kluge avoids shipping the cscript with a parfile.  
91 - """  
92 -  
93 - # Set parfile name  
94 - parfile = self.name+".par"  
95 -  
96 - try:  
97 - pars = gammalib.GApplicationPars(parfile)  
98 - except:  
99 - # Signal that parfile was not found  
100 - print("Parfile \""+parfile+"\" not found. Create default parfile.")  
101 -  
102 - # Check for folders and possible/available configurations  
103 - chains = os.listdir(self.datapath)  
104 -  
105 - # Initialise folder dictionary  
106 - folders = {}  
107 -  
108 - # Loop over possible analysis chains  
109 - for chain in chains:  
110 -  
111 - # Omit hidden files/folders  
112 - if chain[0]=="." or not os.path.isdir(self.datapath+"/"+chain):  
113 - continue  
114 -  
115 - # Initialise dictionary per chain  
116 - folders[chain] = {}  
117 -  
118 - # Read available DST versions  
119 - dstversions = os.listdir(self.datapath+"/"+chain)  
120 -  
121 - # Loop over DST versions  
122 - for dstver in dstversions:  
123 -  
124 - # Omit hidden files/folders  
125 - if dstver[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver):  
126 - continue  
127 -  
128 - # Initialise available cut configurations  
129 - folders[chain][dstver] = []  
130 - configs = os.listdir(self.datapath+"/"+chain+"/"+dstver)  
131 -  
132 - # Loop over configurations  
133 - for config in configs:  
134 -  
135 - # Omit hidden files/folders  
136 - if config[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver+"/"+config):  
137 - continue  
138 -  
139 - # Append configuration to corresponding dictionary  
140 - folders[chain][dstver].append(config)  
141 -  
142 - # Check for appropriate folder structure  
143 - checked_folders = self.check_structure(folders)  
144 - if not len(checked_folders):  
145 - raise KeyError("Appropriate folder structure not found in $HESSFITS. The correct structure can be found on HESS internal documentation")  
146 -  
147 - # Create default parfile  
148 - pars = gammalib.GApplicationPars()  
149 -  
150 - # Create parameters for each available config  
151 - avchains = ""  
152 - for chain in checked_folders:  
153 - avchains += chain +"|"  
154 - avdsts = ""  
155 - for dstver in checked_folders[chain]:  
156 - avdsts+=dstver+"|"  
157 - avconfigs = ""  
158 - for config in checked_folders[chain][dstver]:  
159 - avconfigs += config+"|"  
160 - pars.append(gammalib.GApplicationPar(chain+"_"+dstver+"_config","s","a",checked_folders[chain][dstver][0],avconfigs[:-1],"","Analysis configuration"))  
161 - pars.append(gammalib.GApplicationPar(chain+"_dstver","s","a",checked_folders[chain].keys()[0],avdsts[:-1],"","DST version"))  
162 - pars.append(gammalib.GApplicationPar("chain","s","a",checked_folders.keys()[0],avchains[:-1],"","Analysis chain"))  
163 - pars.append(gammalib.GApplicationPar("runlistfile","f","a","runlist.lis","","","Runlist file"))  
164 - pars.append(gammalib.GApplicationPar("inmodel","f","h","NONE","","","Input model XML file (optional)"))  
165 - pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile"))  
166 - pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file"))  
167 - pars.append(gammalib.GApplicationPar("datapath","s","h",self.datapath,"","","Data directoy"))  
168 - pars.append(gammalib.GApplicationPar("bkgpars","i","a","1","","","Number of free parameters per background model"))  
169 - pars.append(gammalib.GApplicationPar("usetrig","b","h","yes","yes|no","","Pre-estimate background normalisation? Recommended for binned analysis"))  
170 - pars.append(gammalib.GApplicationPar("bkgtype","s","a","irf","gauss|irf|aeff","","What background model should be used"))  
171 - pars.append_standard()  
172 - pars.save(parfile)  
173 -  
174 - # Return  
175 - return  
176 -  
177 - def get_parameters(self):  
178 - """  
179 - Get parameters from parfile and setup the observation.  
180 - """  
181 - # Get parameters  
182 - self.runlistfile = self["runlistfile"].filename()  
183 -  
184 - # Read models if user input  
185 - self.inmodel = self["inmodel"].filename()  
186 - if not (self.inmodel == "NONE" or self.inmodel == ""):  
187 - self.inmodels = gammalib.GModels(self.inmodel)  
188 -  
189 - # Observation outfile  
190 - self.outobs = self["outobs"].filename()  
191 -  
192 - # Output model file  
193 - self.outmodel = self["outmodel"].filename()  
194 -  
195 - # Read hidden parameters  
196 - self.bkgpars = self["bkgpars"].integer()  
197 - self.bkgtype = self["bkgtype"].string()  
198 - self.usetrig = self["usetrig"].boolean()  
199 - self.datapath = self["datapath"].string()  
200 -  
201 - # Read run configuration parameters  
202 - self.chain = self["chain"].string()  
203 - self.dstver = self[self.chain+"_dstver"].string()  
204 - self.config = self[self.chain+"_"+self.dstver+"_config"].string()  
205 -  
206 - # Build directory name prefix  
207 - self.directory = self.datapath+"/"  
208 - self.directory += self.chain +"/"  
209 - self.directory += self.dstver +"/"  
210 - self.directory += self.config +"/"  
211 -  
212 - if not os.path.isdir(self.directory):  
213 - raise gammalib.Exception.runtime_error("Folder structure not found. Please check your HESSFITS environment variable")  
214 -  
215 - self.m_log = False # Logging in client tools  
216 - self.m_debug = False # Debugging in client tools  
217 - self.m_clobber = self["clobber"].boolean()  
218 -  
219 - # Return  
220 - return  
221 -  
222 - def check_structure(self,folders):  
223 -  
224 - # Initialise valid folders  
225 - newfolders = {}  
226 -  
227 - # Check if we have found possible chains  
228 - if len(folders):  
229 -  
230 - # Loop over chains  
231 - for chain in folders:  
232 -  
233 - # Loop over dst versions  
234 - for dst in folders[chain]:  
235 -  
236 - # Loop over configs:  
237 - for config in folders[chain][dst]:  
238 -  
239 - # check if config is valid  
240 - path = self.datapath+"/"+chain+"/"+dst+"/"+config  
241 - if os.path.isdir(path) and len(glob.glob(path+"/run*")):  
242 - if not chain in newfolders:  
243 - newfolders[chain] = {}  
244 - if not dst in newfolders[chain]:  
245 - newfolders[chain][dst] = []  
246 - newfolders[chain][dst].append(config)  
247 -  
248 - # Return valid configs  
249 - return newfolders  
250 -  
251 - def background_spectrum(self, run, prefactor, index, aefffile = ""):  
252 -  
253 - scale = pow(10,math.floor(math.log10(abs(prefactor))))  
254 -  
255 - if index == 0.0 and self.bkgpars <= 1:  
256 - spec = gammalib.GModelSpectralConst()  
257 - spec["Value"].min(0.1*scale)  
258 - spec["Value"].max(10*scale)  
259 - spec["Value"].value(prefactor)  
260 - spec["Value"].scale(scale)  
261 - if self.bkgpars == 0:  
262 - spec["Value"].fix()  
263 - else:  
264 - spec["Value"].free()  
265 -  
266 - else:  
267 -  
268 - # Create power law if number of free parameters are less 2 or less  
269 - if self.bkgpars <= 2:  
270 - e = gammalib.GEnergy(1.0,"TeV")  
271 - spec = gammalib.GModelSpectralPlaw(prefactor, index, e)  
272 -  
273 - # Set parameter ranges  
274 - spec[0].scale(scale)  
275 - spec[0].min(0.01 * scale)  
276 - spec[0].max(100.0 * scale)  
277 - spec[1].scale(1)  
278 - spec[1].min(-5.0)  
279 - spec[1].max(5.0)  
280 -  
281 - # Set number of free parameters  
282 - if self.bkgpars == 0:  
283 - spec[0].fix()  
284 - spec[1].fix()  
285 - elif self.bkgpars == 1:  
286 - spec[0].free()  
287 - spec[1].fix()  
288 - else:  
289 - spec[0].free()  
290 - spec[1].free()  
291 -  
292 - else:  
293 -  
294 - # Use several energy-dependent scaling factors  
295 - fits = gammalib.GFits(aefffile)  
296 - emin = fits["EFFECTIVE AREA"].real("LO_THRES")  
297 - emax = fits["EFFECTIVE AREA"].real("HI_THRES")  
298 - fits.close()  
299 -  
300 - # Create reference powerlaw  
301 - plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,"TeV"))  
302 -  
303 - # Create spectral model and energy values  
304 - spec = gammalib.GModelSpectralNodes()  
305 - bounds = gammalib.GEbounds(self.bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)  
306 - for i in range(bounds.size()):  
307 - energy = bounds.elogmean(i)  
308 - spec.append(energy, plaw.eval(energy, gammalib.GTime()))  
309 - for par in spec:  
310 -  
311 - if "Energy" in par.name():  
312 - par.fix()  
313 - elif "Intensity" in par.name():  
314 - parscale = pow(10,math.floor(math.log10(abs(par.value()))))  
315 - par.scale(parscale)  
316 - par.min(0.01*parscale)  
317 - par.max(100.0*parscale)  
318 -  
319 -  
320 - return spec  
321 -  
322 - def hess_background(self, run, aefffile, trgrate, zenith):  
323 -  
324 - # Set IrfBackground  
325 - if self.bkgtype == "irf":  
326 - prefactor = 1.0  
327 - index = 0.0  
328 - if self.usetrig and trgrate > 0.0:  
329 - # Use value from lookup  
330 - prefactor = 0.00427233352611 * trgrate + 0.181244590076  
331 - spec = self.background_spectrum(run, prefactor, index, aefffile)  
332 - # Create background model instance  
333 - bck = gammalib.GCTAModelIrfBackground(spec)  
334 -  
335 - # Set AeffBackground  
336 - elif self.bkgtype == "aeff":  
337 - prefactor = 5e-14  
338 - index = -2.5  
339 - # Use values from lookup if possible  
340 - if self.usetrig and trgrate > 0.0:  
341 - prefactor = 1.52976410372e-16 * trgrate + 2.71805547108e-14  
342 - if zenith > 0.0:  
343 - index = 0.000271928721732 * zenith * zenith -0.0029602286453 * zenith -2.76404600416  
344 -  
345 - spec = self.background_spectrum(run, prefactor, index, aefffile)  
346 -  
347 - # Create background model instance  
348 - bck = gammalib.GCTAModelAeffBackground(spec)  
349 -  
350 - # Set Gaussian Background  
351 - elif self.bkgtype == "gauss":  
352 - spec = self.background_spectrum(run, 1e-4, -1.8, aefffile)  
353 - radial = gammalib.GCTAModelRadialGauss(2.5)  
354 - bck = gammalib.GCTAModelRadialAcceptance(radial, spec)  
355 -  
356 - else:  
357 - sys.exit("Background type \""+self.bkgtype+"\" unsupported")  
358 -  
359 - # Copy model  
360 - model = bck.clone()  
361 -  
362 - # Assign specific run id  
363 - model.ids(str(int(run)))  
364 -  
365 - # Assign instrument  
366 - model.instruments("HESS")  
367 -  
368 - # Set name (arbitrary)  
369 - model.name("bkg_"+run)  
370 -  
371 - # Turn off TS calculation for background model  
372 - model.tscalc(False)  
373 -  
374 - # Return model  
375 - return model  
376 -  
377 -  
378 - def erange(self):  
379 - """  
380 - Returns runlist energy range  
381 - """  
382 - return self.m_erange  
383 -  
384 - def execute(self):  
385 - """  
386 - Execute the script.  
387 - """  
388 - # Run the script  
389 - self.run()  
390 -  
391 - # Save residual map  
392 - self.save(self.outobs, self.m_clobber)  
393 -  
394 - # Return  
395 - return  
396 -  
397 - # Get runfolder from runnumber (subject to defined folder structure)  
398 - def get_runfolder(self,runnr):  
399 -  
400 - start = runnr[:-2]+"00"  
401 - if not int(runnr[-3]) % 2 == 0:  
402 - start = str(int(start)-100)  
403 - if len(start) == 5:  
404 - start = "0"+start  
405 - int_end = int(start)+199  
406 - end = str(int_end)  
407 - if len(end) == 5:  
408 - end = "0"+ end  
409 -  
410 - runfold = runnr  
411 - if len(runfold) == 5:  
412 - runfold = "0"+runfold  
413 -  
414 - return self.directory+"/run"+start+"-"+end+"/run"+runfold+"/"  
415 -  
416 - def run(self):  
417 - """  
418 - Run the script.  
419 - """  
420 - # Switch screen logging on in debug mode  
421 - if self.logDebug():  
422 - self.log.cout(True)  
423 -  
424 - # Get parameters  
425 - self.get_parameters()  
426 -  
427 - # Write input parameters into logger  
428 - if self.logTerse():  
429 - self.log_parameters()  
430 - self.log("\n")  
431 -  
432 - # Initialise empty models  
433 - self.models = gammalib.GModels()  
434 -  
435 - # Initialise empty xml file and append an observation list  
436 - self.xml = gammalib.GXml()  
437 - self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\""))  
438 - lib = self.xml.element("observation_list", 0)  
439 -  
440 - # Read runlist from file  
441 - runlist = []  
442 - runfile = open(self.runlistfile)  
443 - for line in runfile.readlines():  
444 - if len(line) == 0:  
445 - continue  
446 - if line[0] == '#':  
447 - continue  
448 - if len(line.split()) > 0:  
449 - runlist.append(line.split()[0])  
450 - else:  
451 - runlist.append('')  
452 - runfile.close()  
453 -  
454 - # Loop over runs  
455 - # Log output  
456 - if self.logTerse():  
457 - self.log("\n")  
458 - self.log.header1("Looping over runs")  
459 -  
460 - # Initialise energy range values for logging  
461 - runlist_emin = 100.0  
462 - runlist_emax = 0.0  
463 -  
464 - for runnr in runlist:  
465 -  
466 - # Unify run number string length  
467 - run = runnr  
468 - if len(run) == 5:  
469 - run = "0"+run  
470 -  
471 - # Get folder in which run files should be located  
472 - runfolder = self.get_runfolder(run)  
473 -  
474 - # construct filenames  
475 - eventfile = runfolder + "hess_events_"+run+".fits.gz"  
476 - aefffile = runfolder + "hess_aeff_2d_"+run+".fits.gz"  
477 - psffile = runfolder + "hess_psf_king_"+run+".fits.gz"  
478 - edispfile = runfolder + "hess_edisp_2d_"+run+".fits.gz"  
479 - bgfile = runfolder + "hess_bkg_offruns_"+run+".fits.gz"  
480 -  
481 - # Check for existence of files  
482 - msg = ""  
483 - skip = False  
484 - if not os.path.isdir(runfolder):  
485 - msg = "Run "+str(int(run))+" does not exist in FITS format"  
486 - skip = True  
487 - elif not os.path.isfile(eventfile):  
488 - msg = "Run "+str(int(run))+" has no eventfile - Run is skipped"  
489 - skip = True  
490 - elif not os.path.isfile(aefffile):  
491 - msg = "Run "+str(int(run))+" has no effective area - Run is skipped"  
492 - skip = True  
493 - elif not os.path.isfile(psffile):  
494 - msg = "Run "+str(int(run))+" has no PSF - Run is skipped"  
495 - skip = True  
496 - elif not os.path.isfile(edispfile):  
497 - msg = "Warning: Run "+str(int(run))+" has no energy dispersion - usage of energy dispersion not possible"  
498 - skip = False  
499 - edispfile = ""  
500 - elif not os.path.isfile(bgfile):  
501 - msg = "Warning: Run "+str(int(run))+" has no background file - use of IRF background not possible"  
502 - skip = False  
503 - bgfile = ""  
504 - else:  
505 - msg = "Adding run "+str(int(run))  
506 -  
507 - # Log ouput and skip run if not available  
508 - if self.logTerse():  
509 - self.log(msg+"\n")  
510 - if skip:  
511 - continue  
512 - if self.logExplicit():  
513 - self.log(" Event file: "+eventfile+"\n")  
514 - self.log(" Effective area: "+aefffile+"\n")  
515 - self.log(" Point spread function: "+psffile+"\n")  
516 - self.log(" Energy dispersion: "+edispfile+"\n")  
517 - self.log(" Background file: "+bgfile+"\n")  
518 -  
519 - # Open fits file to determine the observation name  
520 - fits = gammalib.GFits(eventfile)  
521 - object_name = fits["EVENTS"].string("OBJECT")  
522 -  
523 - # Retrieve trigger rate and zenith if available  
524 - trgrate = 0.0  
525 - zenith = 0.0  
526 - if fits["EVENTS"].has_card("ZTRGRATE"):  
527 - trgrate = fits["EVENTS"].real("ZTRGRATE")  
528 - if fits["EVENTS"].has_card("ALT_PNT"):  
529 - zenith = 90.-fits["EVENTS"].real("ALT_PNT")  
530 -  
531 - # Close FITS file  
532 - fits.close()  
533 -  
534 - # Open effective area file to look for threshold  
535 - aeff_fits = gammalib.GFits(aefffile)  
536 - run_emin = aeff_fits["EFFECTIVE AREA"].real("LO_THRES")  
537 - run_emax = aeff_fits["EFFECTIVE AREA"].real("HI_THRES")  
538 - if run_emin < runlist_emin:  
539 - runlist_emin = run_emin  
540 - if run_emax > runlist_emax:  
541 - runlist_emax = run_emax  
542 - aeff_fits.close()  
543 -  
544 - # Append observation to XML and set attributes  
545 - obs = lib.append("observation");  
546 - obs.attribute("name", object_name);  
547 - obs.attribute("id", str(int(run)));  
548 - obs.attribute("instrument", "HESS");  
549 -  
550 - # Append event file  
551 - par = gammalib.GXmlElement("parameter name=\"EventList\"")  
552 - par.attribute("file", eventfile);  
553 - obs.append(par)  
554 -  
555 - # Append effective area  
556 - aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"")  
557 - aeff.attribute("file",aefffile)  
558 -  
559 - # Append PSF  
560 - psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"")  
561 - psf.attribute("file",psffile)  
562 -  
563 - # Append energy dispersion  
564 - edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"")  
565 - edisp.attribute("file",edispfile)  
566 -  
567 - # Append background  
568 - bck = gammalib.GXmlElement("parameter name=\"Background\"")  
569 - bck.attribute("file",bgfile)  
570 -  
571 - # assign events and IRFs to observation  
572 - obs.append(aeff)  
573 - obs.append(psf)  
574 - obs.append(edisp)  
575 - obs.append(bck)  
576 -  
577 - # Append instrumental background model  
578 - self.models.append(self.hess_background(run, aefffile, trgrate, zenith))  
579 -  
580 - # Continue only if there are observations available  
581 - if lib.size():  
582 -  
583 - # Store energy range as member  
584 - self.m_erange = gammalib.GEbounds(gammalib.GEnergy(runlist_emin,"TeV"),gammalib.GEnergy(runlist_emax,"TeV"))  
585 -  
586 - # Write energy range into log file  
587 - if self.logTerse():  
588 - self.log("\n")  
589 - self.log("Eenrgy range of obervation list: ")  
590 - self.log(("% 3.2f" % runlist_emin)+" - "+("% 3.2f" % runlist_emax)+" TeV")  
591 -  
592 - else:  
593 - self.log.header2("WARNING: No observation from given runlist available")  
594 -  
595 - # Append models provided by 'inmodels' if necessary  
596 - if not self.inmodels == None:  
597 - if self.logTerse():  
598 - self.log("\n")  
599 - self.log.header1("Appending models")  
600 -  
601 - # Loop over input models  
602 - for model in self.inmodels:  
603 - if self.logTerse():  
604 - self.log.header3("Adding model \""+model.name()+"\"")  
605 - self.models.append(model)  
606 -  
607 - # Return  
608 - return  
609 -  
610 - def save(self,outfile,clobber):  
611 -  
612 - # Save observation XML file  
613 - self.xml.save(outfile)  
614 -  
615 - # Save model XML file  
616 - self.models.save(self.outmodel)  
617 -  
618 -  
619 -# ======================== #  
620 -# Main routine entry point #  
621 -# ======================== #  
622 -if __name__ == '__main__':  
623 - """  
624 - Generates observation file.  
625 - """  
626 - # Create instance of application  
627 - app = cshessobs(sys.argv)  
628 -  
629 - # Open logfile  
630 - app.logFileOpen()  
631 -  
632 - # Execute application  
633 - app.execute()  
cscripts/csiactdata.py 0 → 100755
  View file @ea3e9e1
  1 +#!/usr/bin/env python
  2 +# ==========================================================================
  3 +# Inspection of IACT data storage
  4 +#
  5 +# Copyright (C) 2015 Michael Mayer
  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 +import gammalib
  22 +import ctools
  23 +import sys
  24 +import os
  25 +import json
  26 +
  27 +
  28 +# =============== #
  29 +# cshessobs class #
  30 +# =============== #
  31 +class csiactdata(ctools.cscript):
  32 + """
  33 + This script inspects the available FITS data storage on the user machine and writes
  34 + information about the available FITS configurations into a log file or prints it
  35 + on the screen. These information can be used as input for 'csiactobs'.
  36 + """
  37 + def __init__(self, *argv):
  38 + """
  39 + Constructor.
  40 + """
  41 + # Set name
  42 + self.name = "csiactdata"
  43 + self.version = "1.0.0"
  44 + self.datapath = os.getenv("VHEFITS","")
  45 +
  46 + # Make sure that parfile exists
  47 + file = 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)
  60 +
  61 + # Return
  62 + return
  63 +
  64 + def __del__(self):
  65 + """
  66 + Destructor.
  67 + """
  68 + # Return
  69 + return
  70 +
  71 + def parfile(self):
  72 + """
  73 + Check if parfile exists. If parfile does not exist then create a
  74 + default parfile. This kluge avoids shipping the cscript with a parfile.
  75 + """
  76 + # Set parfile name
  77 + parfile = self.name+".par"
  78 +
  79 + try:
  80 + pars = gammalib.GApplicationPars(parfile)
  81 + except:
  82 + # Signal if parfile was not found
  83 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  84 +
  85 + # Create default parfile
  86 + pars = gammalib.GApplicationPars()
  87 + pars.append(gammalib.GApplicationPar("datapath","s","a",self.datapath,"","","Path were data is located"))
  88 + pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file"))
  89 + pars.append_standard()
  90 + pars.save(parfile)
  91 +
  92 + # Return
  93 + return
  94 +
  95 + def get_parameters(self):
  96 + """
  97 + Get parameters from parfile and setup the observation.
  98 + """
  99 +
  100 + # Get datapath if not already set
  101 + if self.datapath == "":
  102 + self.datapath = self["datapath"].string()
  103 +
  104 + # Expand environment
  105 + self.datapath = gammalib.expand_env(self.datapath)
  106 +
  107 + # Get filename of master index file
  108 + self.m_master_indx = self["master_indx"].string()
  109 + self.m_master_file = os.path.join(self.datapath, self.m_master_indx)
  110 +
  111 + # Return
  112 + return
  113 +
  114 +
  115 + def execute(self):
  116 + """
  117 + Execute the script.
  118 + """
  119 + # Run the script
  120 + self.run()
  121 +
  122 + # Return
  123 + return
  124 +
  125 + def run(self):
  126 + """
  127 + Run the script.
  128 + """
  129 + # Switch screen logging on
  130 + self.log.cout(True)
  131 +
  132 + # Get parameters
  133 + self.get_parameters()
  134 +
  135 + # Write input parameters into logger
  136 + if self.logTerse():
  137 + self.log_parameters()
  138 + self.log("\n")
  139 +
  140 + self.log("\n")
  141 + self.log.header1("Data storage entries")
  142 + self.log("\n")
  143 + self.log.parformat("Master index file")
  144 + self.log(self.m_master_file)
  145 + self.log("\n")
  146 + self.log("\n")
  147 +
  148 +
  149 + json_data = open(self.m_master_file).read()
  150 + data = json.loads(json_data)
  151 + configs = data["datasets"]
  152 +
  153 + # Loop over configs and display unavailable storage first
  154 + for config in configs:
  155 +
  156 + # Create hdu and obs index files
  157 + hdu = os.path.join(self.datapath, config["hduindx"])
  158 + obs = os.path.join(self.datapath, config["obsindx"])
  159 +
  160 + # Check if index files are available
  161 + if not (gammalib.is_fits(str(hdu)) and gammalib.is_fits(str(obs))):
  162 + self.log.parformat(str(config["name"]))
  163 + self.log("Not available")
  164 + self.log("\n")
  165 +
  166 + # Loop over configs and log available configs
  167 + for config in configs:
  168 +
  169 + # Create hdu and obs index files
  170 + hdu = os.path.join(self.datapath, config["hduindx"])
  171 + obs = os.path.join(self.datapath, config["obsindx"])
  172 +
  173 + # Check if index files are available
  174 + if gammalib.is_fits(str(hdu)) and gammalib.is_fits(str(obs)):
  175 +
  176 + # Log data information
  177 + self.log("\n")
  178 + self.log.header3(str(config["name"]))
  179 + for key in config:
  180 + self.log.parformat(str(key))
  181 + self.log(str(config[key]))
  182 + self.log("\n")
  183 +
  184 + # Return
  185 + return
  186 +
  187 +
  188 +# ======================== #
  189 +# Main routine entry point #
  190 +# ======================== #
  191 +if __name__ == '__main__':
  192 + """
  193 + Shows data storage information
  194 + """
  195 + # Create instance of application
  196 + app = csiactdata(sys.argv)
  197 +
  198 + # Open logfile
  199 + app.logFileOpen()
  200 +
  201 + # Execute application
  202 + app.execute()
cscripts/csiactobs.py 0 → 100755
  View file @ea3e9e1
  1 +#!/usr/bin/env python
  2 +# ==========================================================================
  3 +# Generation of an IACT observation definition file.
  4 +#
  5 +# Copyright (C) 2015 Michael Mayer
  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 +import gammalib
  22 +import ctools
  23 +import sys
  24 +import os
  25 +import json
  26 +
  27 +
  28 +# =============== #
  29 +# cshessobs class #
  30 +# =============== #
  31 +class csiactobs(ctools.cscript):
  32 + """
  33 + This class implements the creation of a observation xml file for IACT
  34 + data analysis. This class is dedicated for use inside a IACT Collaboration,
  35 + i.e. it can only be used if you have access to IACT data in FITS format.
  36 + The FITS data has to be structured and in the format described here:
  37 + http://gamma-astro-data-formats.readthedocs.org/en/latest/
  38 + """
  39 + def __init__(self, *argv):
  40 + """
  41 + Constructor.
  42 + """
  43 + # Set name
  44 + self.name = "csiactobs"
  45 + self.version = "1.0.0"
  46 + self.m_erange = gammalib.GEbounds()
  47 +
  48 + self.datapath = os.getenv("VHEFITS","")
  49 + print argv
  50 + # Initialise some members
  51 + self.inmodels = None
  52 + self.outobs = "obs.xml"
  53 +
  54 + # Make sure that parfile exists
  55 + file = self.parfile()
  56 +
  57 + # Initialise application
  58 + if len(argv) == 0:
  59 + ctools.cscript.__init__(self, self.name, self.version)
  60 + elif len(argv) == 1:
  61 + ctools.cscript.__init__(self, self.name, self.version, *argv)
  62 + else:
  63 + raise TypeError("Invalid number of arguments given.")
  64 +
  65 + # Set logger properties
  66 + self.log_header()
  67 + self.log.date(True)
  68 +
  69 + # Return
  70 + return
  71 +
  72 + def __del__(self):
  73 + """
  74 + Destructor.
  75 + """
  76 + # Return
  77 + return
  78 +
  79 + def parfile(self):
  80 + """
  81 + Check if parfile exists. If parfile does not exist then create a
  82 + default parfile. This kluge avoids shipping the cscript with a parfile.
  83 + """
  84 + # Set parfile name
  85 + parfile = self.name+".par"
  86 +
  87 + try:
  88 + pars = gammalib.GApplicationPars(parfile)
  89 + except:
  90 + # Signal if parfile was not found
  91 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  92 +
  93 + # Create default parfile
  94 + pars = gammalib.GApplicationPars()
  95 + pars.append(gammalib.GApplicationPar("datapath","s","a",self.datapath,"","","Path were data is located"))
  96 + pars.append(gammalib.GApplicationPar("prodname","s","a","","","","Data storage name (Run csdsinfo to view your options)"))
  97 + pars.append(gammalib.GApplicationPar("runlistfile","f","a","runlist.lis","","","Runlist file"))
  98 + pars.append(gammalib.GApplicationPar("bkgpars","i","a","1","","","Number of free parameters per background model"))
  99 + pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile"))
  100 + pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file"))
  101 + pars.append(gammalib.GApplicationPar("inmodel","f","h","NONE","","","Input model XML file (optional)"))
  102 + pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file"))
  103 + pars.append(gammalib.GApplicationPar("ev_hiera","s","h","events","","","Hierarchy of event formats"))
  104 + pars.append(gammalib.GApplicationPar("aeff_hiera","s","h","aeff_2d","","","Hierarchy of effective area formats"))
  105 + pars.append(gammalib.GApplicationPar("psf_hiera","s","h","psf_king|psf_3gauss","","","Hierarchy of psf formats"))
  106 + pars.append(gammalib.GApplicationPar("edisp_hiera","s","h","edisp_2d","","","Hierarchy of energy dispersion formats"))
  107 + pars.append(gammalib.GApplicationPar("bkg_hiera","s","h","bkg_3d","","","Hierarchy of background formats"))
  108 + pars.append(gammalib.GApplicationPar("bkg_mod_hiera","s","h","irf","irf|aeff|gauss","","Hierarchy of background models"))
  109 + pars.append(gammalib.GApplicationPar("bkg_gauss_norm","r","h","1e-8","","","Input normalisation for Gaussian background"))
  110 + pars.append(gammalib.GApplicationPar("bkg_gauss_index","r","h","-2.0","","","Input spectral index for Gaussian background"))
  111 + pars.append(gammalib.GApplicationPar("bkg_gauss_sigma","r","h","2.5","","","Input sigma for Gaussian background"))
  112 + pars.append(gammalib.GApplicationPar("bkg_aeff_norm","r","h","1e-14","","","Input normalisation for effective area background"))
  113 + pars.append(gammalib.GApplicationPar("bkg_aeff_index","r","h","-2.0","","","Input spectral index for effective area background"))
  114 + pars.append(gammalib.GApplicationPar("bkg_range_factor","r","h","100.0","","","Factor to determine range of background normalisation"))
  115 +
  116 + pars.append_standard()
  117 + pars.save(parfile)
  118 +
  119 + # Return
  120 + return
  121 +
  122 + def get_parameters(self):
  123 + """
  124 + Get parameters from parfile and setup the observation.
  125 + """
  126 +
  127 + if self.datapath == "":
  128 + self.datapath = self["datapath"].string()
  129 +
  130 + # Expand environment
  131 + self.datapath = gammalib.expand_env(self.datapath)
  132 +
  133 + # Read user parameters
  134 + self.m_prodname = self["prodname"].string()
  135 + self.m_runlistfile = gammalib.expand_env(self["runlistfile"].filename())
  136 + self.m_bkgpars = self["bkgpars"].integer()
  137 +
  138 + # Output model file
  139 + self.outmodel = self["outmodel"].filename()
  140 +
  141 + # Observation outfile
  142 + self.outobs = self["outobs"].filename()
  143 +
  144 + # Master index file name
  145 + self.m_master_indx = self["master_indx"].string()
  146 +
  147 + # Read hierarchy of file loading
  148 + self.m_ev_hiera = self["ev_hiera"].string().split("|")
  149 + self.m_aeff_hiera = self["aeff_hiera"].string().split("|")
  150 + self.m_psf_hiera = self["psf_hiera"].string().split("|")
  151 + self.m_bkg_hiera = self["bkg_hiera"].string().split("|")
  152 + self.m_edisp_hiera = self["edisp_hiera"].string().split("|")
  153 + self.m_bkg_mod_hiera = self["bkg_mod_hiera"].string().split("|")
  154 +
  155 + # Read hidden background parameters
  156 + self.m_bkg_gauss_norm = self["bkg_gauss_norm"].real()
  157 + self.m_bkg_gauss_index = self["bkg_gauss_index"].real()
  158 + self.m_bkg_gauss_sigma = self["bkg_gauss_sigma"].real()
  159 + self.m_bkg_aeff_index = self["bkg_aeff_index"].real()
  160 + self.m_bkg_aeff_norm = self["bkg_aeff_norm"].real()
  161 + self.m_bkg_range_factor = self["bkg_range_factor"].real()
  162 +
  163 + # Open master index file and look for prodname
  164 + master_file = os.path.join(self.datapath, self.m_master_indx)
  165 + if not os.path.isfile(master_file):
  166 + raise RuntimeError("FITS data store not available. No master index file found. Make sure the file is copied from the server and your datapath is set correctly.")
  167 +
  168 + json_data = open(master_file).read()
  169 + data = json.loads(json_data)
  170 +
  171 + if not "datasets" in data:
  172 + raise RuntimeError("Key \"datasets\" not available in master index file.")
  173 +
  174 + configs = data["datasets"]
  175 +
  176 + self.m_hdu_index = self.m_obs_index = ""
  177 +
  178 + for config in configs:
  179 + if self.m_prodname == config["name"]:
  180 + self.m_hdu_index = str(os.path.join(self.datapath, config["hduindx"]))
  181 + self.m_obs_index = str(os.path.join(self.datapath, config["obsindx"]))
  182 + break
  183 +
  184 + if self.m_hdu_index == "" or self.m_obs_index == "":
  185 + raise RuntimeError("FITS data store \""+self.m_prodname+"\" not available. Run csdsinfo to get a list of available storage names")
  186 +
  187 + if not gammalib.is_fits(self.m_hdu_index+"[HDU_INDEX]"):
  188 + raise RuntimeError("HDU index file \""+self.m_hdu_index+"[HDU_INDEX]\" for FITS data store \""+self.m_prodname+"\" not available. Check your master index file or run csdsinfo to get a list of available storage names.")
  189 +
  190 + if gammalib.is_fits(self.m_obs_index+"[OBS_INDEX]"):
  191 + fits = gammalib.GFits(self.m_obs_index)
  192 + self.has_bkg_scale = fits["OBS_INDEX"].contains("BKG_SCALE")
  193 + fits.close()
  194 + else:
  195 + self.has_bkg_scale = False
  196 +
  197 + # Create base data directory from hdu index file location
  198 + self.subdir = os.path.dirname(self.m_hdu_index)
  199 +
  200 + self.m_log = False # Logging in client tools
  201 + self.m_debug = False # Debugging in client tools
  202 + self.m_clobber = self["clobber"].boolean()
  203 +
  204 + # Return
  205 + return
  206 +
  207 + def background_spectrum(self, run, prefactor, index, emin, emax):
  208 +
  209 + # Handle constant spectral model
  210 + if index == 0.0 and self.m_bkgpars <= 1:
  211 + spec = gammalib.GModelSpectralConst()
  212 + spec["Value"].min(prefactor / self.m_bkg_range_factor)
  213 + spec["Value"].max(prefactor * self.m_bkg_range_factor)
  214 + spec["Value"].value(prefactor)
  215 + if self.m_bkgpars == 0:
  216 + spec["Value"].fix()
  217 + else:
  218 + spec["Value"].free()
  219 +
  220 + else:
  221 +
  222 + # Create power law model
  223 + if self.m_bkgpars <= 2:
  224 + e = gammalib.GEnergy(1.0,"TeV")
  225 + spec = gammalib.GModelSpectralPlaw(prefactor, index, e)
  226 +
  227 + # Set parameter ranges
  228 + spec[0].min(0.01 * prefactor)
  229 + spec[0].max(100.0 * prefactor)
  230 + spec[1].scale(1)
  231 + spec[1].min(-5.0)
  232 + spec[1].max(5.0)
  233 +
  234 + # Set number of free parameters
  235 + if self.m_bkgpars == 0:
  236 + spec[0].fix()
  237 + spec[1].fix()
  238 + elif self.m_bkgpars == 1:
  239 + spec[0].free()
  240 + spec[1].fix()
  241 + else:
  242 + spec[0].free()
  243 + spec[1].free()
  244 +
  245 + else:
  246 +
  247 + # Create reference powerlaw
  248 + plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,"TeV"))
  249 +
  250 + # Create spectral model and energy values
  251 + spec = gammalib.GModelSpectralNodes()
  252 + bounds = gammalib.GEbounds(self.bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)
  253 + for i in range(bounds.size()):
  254 + energy = bounds.elogmean(i)
  255 + value = plaw.eval(energy, gammalib.GTime())
  256 + spec.append(energy, value)
  257 + for par in spec:
  258 +
  259 + if "Energy" in par.name():
  260 + par.fix()
  261 + elif "Intensity" in par.name():
  262 + par.min(value / self.m_bkg_range_factor)
  263 + par.max(value * self.m_bkg_range_factor)
  264 +
  265 + # Return spectrum
  266 + return spec
  267 +
  268 + def iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.1, emax=100):
  269 +
  270 + # handle IrfBackground
  271 + if bkgtype == "irf":
  272 +
  273 + prefactor = 1.0
  274 + index = 0.0
  275 + prefactor *= bkg_scale
  276 + spec = self.background_spectrum(obs_id, prefactor, index, emin, emax)
  277 +
  278 + # Create background model instance
  279 + bck = gammalib.GCTAModelIrfBackground(spec)
  280 +
  281 + # Set AeffBackground
  282 + elif self.bkgtype == "aeff":
  283 +
  284 + prefactor = bkg_scale * self.m_bkg_aeff_norm
  285 + spec = self.background_spectrum(obs_id, prefactor, self.m_bkg_aeff_index, emin, emax)
  286 +
  287 + # Create background model instance
  288 + bck = gammalib.GCTAModelAeffBackground(spec)
  289 +
  290 + # Set Gaussian Background
  291 + elif self.bkgtype == "gauss":
  292 +
  293 + prefactor = bkg_scale * self.m_bkg_gauss_norm
  294 + spec = self.background_spectrum(obs_id, prefactor, self.m_bkg_gauss_index, emin, emax)
  295 + radial = gammalib.GCTAModelRadialGauss(self.m_bkg_gauss_sigma)
  296 + bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
  297 +
  298 + else:
  299 + sys.exit("Background type \""+self.bkgtype+"\" unsupported")
  300 +
  301 + # Copy model
  302 + model = bck.clone()
  303 +
  304 + # Assign specific run id
  305 + model.ids(str(obs_id))
  306 +
  307 + # Assign instrument
  308 + model.instruments(telescope)
  309 +
  310 + # Set name (arbitrary)
  311 + model.name("bkg_"+str(obs_id))
  312 +
  313 + # Turn off TS calculation for background model
  314 + model.tscalc(False)
  315 +
  316 + # Return model
  317 + return model
  318 +
  319 +
  320 + def erange(self):
  321 + """
  322 + Returns runlist energy range
  323 + """
  324 + return self.m_erange
  325 +
  326 + def execute(self):
  327 + """
  328 + Execute the script.
  329 + """
  330 + # Run the script
  331 + self.run()
  332 +
  333 + # Save residual map
  334 + self.save(self.outobs, self.m_clobber)
  335 +
  336 + # Return
  337 + return
  338 +
  339 + def run(self):
  340 + """
  341 + Run the script.
  342 + """
  343 + # Switch screen logging on in debug mode
  344 + if self.logDebug():
  345 + self.log.cout(True)
  346 +
  347 + # Get parameters
  348 + self.get_parameters()
  349 +
  350 + # Write input parameters into logger
  351 + if self.logTerse():
  352 + self.log_parameters()
  353 + self.log("\n")
  354 +
  355 + # Initialise empty models
  356 + self.models = gammalib.GModels()
  357 +
  358 + # Initialise empty xml file and append an observation list
  359 + self.xml = gammalib.GXml()
  360 + self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\""))
  361 + lib = self.xml.element("observation_list", 0)
  362 +
  363 + # Read runlist from file
  364 + runlist = []
  365 + runfile = open(self.m_runlistfile)
  366 + for line in runfile.readlines():
  367 + if len(line) == 0:
  368 + continue
  369 + if line[0] == '#':
  370 + continue
  371 + if len(line.split()) > 0:
  372 + runlist.append(line.split()[0])
  373 + else:
  374 + runlist.append('')
  375 + runfile.close()
  376 +
  377 + # Loop over runs
  378 + # Log output
  379 + if self.logTerse():
  380 + self.log("\n")
  381 + self.log.header1("Looping over runs")
  382 +
  383 + # Initialise energy range values for logging
  384 + runlist_emin = 100.0
  385 + runlist_emax = 0.0
  386 +
  387 + for obs_id in runlist:
  388 +
  389 + obs_selection = "[OBS_ID=="+str(obs_id)+"]"
  390 +
  391 + hduindx = gammalib.GFits(self.m_hdu_index+"[HDU_INDEX]"+obs_selection)
  392 + hduindx_hdu = hduindx["HDU_INDEX"]
  393 +
  394 + eventfile = aefffile = psffile = edispfile = bkgfile = ""
  395 + eventhdu = aeffhdu = ""
  396 +
  397 + types = []
  398 + formats = []
  399 + for i in range(hduindx_hdu.nrows()):
  400 + types.append(hduindx_hdu["HDU_TYPE"][i])
  401 + formats.append(hduindx_hdu["HDU_CLASS"][i])
  402 +
  403 + # Handle events
  404 + index = -1
  405 + for version in self.m_ev_hiera:
  406 + n = formats.count(version)
  407 + if n > 0:
  408 + index = formats.index(version)
  409 + break
  410 + if index >= 0:
  411 + eventfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index])
  412 + eventhdu = hduindx_hdu["HDU_NAME"][index]
  413 + eventfile += "["+eventhdu+"]"
  414 + if not gammalib.is_fits(eventfile):
  415 + if self.logTerse():
  416 + self.log("Skipping observation "+str(obs_id)+": eventfile \""+eventfile+"\" not found\n")
  417 + continue
  418 +
  419 + # Handle Aeff
  420 + index = -1
  421 + for version in self.m_aeff_hiera:
  422 + n = formats.count(version)
  423 + if n > 0:
  424 + index = formats.index(version)
  425 + break
  426 + if index >= 0:
  427 + aefffile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index])
  428 + aeffhdu = hduindx_hdu["HDU_NAME"][index]
  429 + aefffile += "["+aeffhdu+"]"
  430 + if not gammalib.is_fits(aefffile):
  431 + if self.logTerse():
  432 + self.log("Skipping observation "+str(obs_id)+": effective area \""+aefffile+"\" not found\n")
  433 + continue
  434 +
  435 + # Handle psf
  436 + index = -1
  437 + for version in self.m_psf_hiera:
  438 + n = formats.count(version)
  439 + if n > 0:
  440 + index = formats.index(version)
  441 + break
  442 + if index >= 0:
  443 + psffile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index])
  444 + psffile += "["+hduindx_hdu["HDU_NAME"][index]+"]"
  445 + if not gammalib.is_fits(psffile):
  446 + if self.logTerse():
  447 + self.log("Skipping observation "+str(obs_id)+": point spread function \""+psffile+"\" not found\n")
  448 + continue
  449 +
  450 + # Handle edisp
  451 + index = -1
  452 + for version in self.m_edisp_hiera:
  453 + n = formats.count(version)
  454 + if n > 0:
  455 + index = formats.index(version)
  456 + break
  457 + if index >= 0:
  458 + edispfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index])
  459 + edispfile += "["+hduindx_hdu["HDU_NAME"][index]+"]"
  460 + if not gammalib.is_fits(edispfile):
  461 + if self.logTerse():
  462 + self.log("Warning: observation "+str(obs_id)+" has no energy dispersion \""+edispfile+"\" information\n")
  463 + edispfile = ""
  464 +
  465 + # Handle background
  466 + index = -1
  467 + bkg_mod_hierarchy = self.m_bkg_mod_hiera
  468 + for version in self.m_bkg_hiera:
  469 + n = formats.count(version)
  470 + if n > 0:
  471 + index = formats.index(version)
  472 + break
  473 + if index >= 0:
  474 + bkgfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index])
  475 + bkgfile += "["+hduindx_hdu["HDU_NAME"][index]+"]"
  476 + if not gammalib.is_fits(bkgfile):
  477 + if self.logTerse():
  478 + self.log("Warning: observation "+str(obs_id)+" has no background information in \""+bkgfile+"\". IRF background cannot be used\n")
  479 + bkgfile = ""
  480 + bkg_mod_hierarchy.remove("irf")
  481 +
  482 + # Handle background scale information if available
  483 + bkg_scale = 1.0
  484 + if self.has_bkg_scale:
  485 + obsindx = gammalib.GFits(self.m_obs_index+"[OBS_INDEX]"+obs_selection)
  486 + bkg_scale = obsindx["OBS_INDEX"]["BKG_SCALE"][0]
  487 + obsindx.close()
  488 +
  489 + # Open fits file to determine the observation name
  490 + fits = gammalib.GFits(eventfile)
  491 + events = fits[eventhdu]
  492 + object_name = events.string("OBJECT")
  493 + telescope = events.string("TELESCOP")
  494 +
  495 + # Close FITS file
  496 + fits.close()
  497 +
  498 + # Open effective area file to look for threshold
  499 + aeff_fits = gammalib.GFits(aefffile)
  500 + aeff_table = aeff_fits[aeffhdu]
  501 +
  502 + # Initialise energy range in case no threshold is available
  503 + run_emin = 0.1
  504 + run_emax = 100.0
  505 +
  506 + # Set energy range from heady keywords
  507 + if aeff_table.has_card("LO_THRES") and aeff_table.has_card["HI_THRES"]:
  508 + run_emin = aeff_table.real("LO_THRES")
  509 + run_emax = aeff_table.real("HI_THRES")
  510 +
  511 + # Close Aeff fits file
  512 + aeff_fits.close()
  513 +
  514 + # Adapt runlist energy range if necessary
  515 + if run_emin < runlist_emin:
  516 + runlist_emin = run_emin
  517 + if run_emax > runlist_emax:
  518 + runlist_emax = run_emax
  519 +
  520 + # Logging
  521 + if self.logTerse():
  522 + self.log("Adding observation "+str(obs_id)+" (\""+object_name+"\")\n")
  523 +
  524 + if self.logExplicit():
  525 + self.log(" Event file: "+eventfile+"\n")
  526 + self.log(" Effective area: "+aefffile+"\n")
  527 + self.log(" Point spread function: "+psffile+"\n")
  528 + self.log(" Energy dispersion: "+edispfile+"\n")
  529 + self.log(" Background file: "+bkgfile+"\n")
  530 + if self.has_bkg_scale:
  531 + self.log(" Background scale: "+str(bkg_scale)+"\n")
  532 + self.log("\n")
  533 +
  534 + # Append observation to XML and set attributes
  535 + obs = lib.append("observation");
  536 + obs.attribute("name", object_name);
  537 + obs.attribute("id", str(obs_id));
  538 + obs.attribute("instrument", telescope);
  539 +
  540 + # Append event file
  541 + ev = gammalib.GXmlElement("parameter name=\"EventList\"")
  542 + ev.attribute("file", eventfile);
  543 +
  544 + # Append effective area
  545 + aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"")
  546 + aeff.attribute("file",aefffile)
  547 +
  548 + # Append PSF
  549 + psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"")
  550 + psf.attribute("file",psffile)
  551 +
  552 + # Append energy dispersion
  553 + edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"")
  554 + edisp.attribute("file",edispfile)
  555 +
  556 + # Append background
  557 + bck = gammalib.GXmlElement("parameter name=\"Background\"")
  558 + bck.attribute("file",bkgfile)
  559 +
  560 + # assign events and IRFs to observation
  561 + obs.append(ev)
  562 + obs.append(aeff)
  563 + obs.append(psf)
  564 + obs.append(edisp)
  565 + obs.append(bck)
  566 +
  567 + # Append instrumental background model
  568 + self.models.append(self.iact_background(telescope, obs_id, bkg_scale, bkg_mod_hierarchy[0], run_emin, run_emax))
  569 +
  570 + # Continue only if there are observations available
  571 + if lib.size():
  572 +
  573 + # Store energy range as member
  574 + self.m_erange = gammalib.GEbounds(gammalib.GEnergy(runlist_emin,"TeV"),gammalib.GEnergy(runlist_emax,"TeV"))
  575 +
  576 + # Write energy range into log file
  577 + if self.logTerse():
  578 + self.log("\n")
  579 + self.log("Eenrgy range of obervation list: ")
  580 + self.log(("% 3.2f" % runlist_emin)+" - "+("% 3.2f" % runlist_emax)+" TeV")
  581 + self.log("\n")
  582 +
  583 + else:
  584 + self.log.header2("WARNING: No observation from given runlist available")
  585 +
  586 + # Append models provided by 'inmodels' if necessary
  587 + if not self.inmodels == None:
  588 + if self.logTerse():
  589 + self.log("\n")
  590 + self.log.header1("Appending models")
  591 +
  592 + # Loop over input models
  593 + for model in self.inmodels:
  594 + if self.logTerse():
  595 + self.log.header3("Adding model \""+model.name()+"\"")
  596 + self.models.append(model)
  597 +
  598 + # Return
  599 + return
  600 +
  601 + def save(self,outfile,clobber):
  602 +
  603 + # Save observation XML file
  604 + self.xml.save(outfile)
  605 +
  606 + # Save model XML file
  607 + self.models.save(self.outmodel)
  608 +
  609 +
  610 +# ======================== #
  611 +# Main routine entry point #
  612 +# ======================== #
  613 +if __name__ == '__main__':
  614 + """
  615 + Generates observation file.
  616 + """
  617 + # Create instance of application
  618 + app = csiactobs(sys.argv)
  619 +
  620 + # Open logfile
  621 + app.logFileOpen()
  622 +
  623 + # Execute application
  624 + app.execute()
test/iactdata/hdu-index.fits.gz 0 → 100644
View file @ea3e9e1
No preview for this file type
test/iactdata/master.json 0 → 100644
  View file @ea3e9e1
  1 +{
  2 + "datasets": [
  3 + {
  4 + "name": "unit-test",
  5 + "obsindx": "obs-index.fits.gz",
  6 + "hduindx": "hdu-index.fits.gz"
  7 + }
  8 + ]
  9 +}
test/iactdata/obs-index.fits.gz 0 → 100644
View file @ea3e9e1
No preview for this file type
test/iactdata/runlist.txt 0 → 100644
  View file @ea3e9e1
  1 +0
test/test_cscripts.sh.in
  View file @ea3e9e1
@@ -114,6 +114,28 @@ else @@ -114,6 +114,28 @@ else
114 fi 114 fi
115 $ECHO " ok" 115 $ECHO " ok"
116 116
  117 +#
  118 +# Test csiactobs
  119 +# ==============
  120 +$ECHO -n "Test csiactobs: "
  121 +chmod +x @abs_top_srcdir@/cscripts/csiactobs.py
  122 +@abs_top_srcdir@/cscripts/csiactobs.py \
  123 + datapath="iactdata/" \
  124 + prodname="unit-test" \
  125 + runlistfile="iactdata/runlist.txt" \
  126 + bkgpars="1" \
  127 + outobs="obs.xml" \
  128 + outmodel="bkgmodels.xml" > csiactobs.log
  129 +$ECHO -n "."
  130 +if [ -s "obs.xml" ]
  131 +then
  132 + $ECHO -n "."
  133 +else
  134 + $ECHO " obs.xml file is not valid"
  135 + exit 1
  136 +fi
  137 +$ECHO " ok"
  138 +
117 139
118 # 140 #
119 # Test csmodelinfo 141 # Test csmodelinfo