Wrijupan Bhattacharyya / ctools

Commit e5a3b78c91ea12018cfb96c62c7f8c716363443d

Authored by Michael Mayer
1 parent 318346b4

added cshessobs script to build observation xml files from a given runlist

Inline Side-by-side
Showing 2 changed files with 518 additions and 2 deletions
scripts/Makefile.am
  View file @e5a3b78
... ... @@ -5,14 +5,14 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \
5 5 tsmerge.py
6 6  
7 7 # Scripts to be distributed in pythondir
8   -pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
  8 +pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py cshessobs.py \
9 9 csobsdef.py csspec.py \
10 10 lsmodel.py \
11 11 obsutils.py
12 12  
13 13 # Scripts for which symbolic links will be created in bindir
14 14 linkscripts = cspull cssens cstsdist cscaldb csresmap \
15   - csobsdef csspec \
  15 + csobsdef csspec cshessobs\
16 16 lsmodel
17 17  
18 18 # Distribute scripts
... ...
scripts/cshessobs.py 0 → 100755
  View file @e5a3b78
  1 +#!/usr/bin/env python
  2 +import gammalib
  3 +import ctools
  4 +import sys
  5 +import os
  6 +import glob
  7 +
  8 +# =============== #
  9 +# cshessobs class #
  10 +# =============== #
  11 +class cshessobs(ctools.cscript):
  12 + """
  13 + This class implements the creation of a observation xml file. It derives
  14 + from the GammaLib::GApplication class which provides support for parameter
  15 + files, command line arguments, and logging. In that way the Python
  16 + script behaves just as a regular ctool. This script is dedicated for use inside
  17 + the HESS Collaboration, i.e. it can only be used if you have access to HESS data
  18 + in FITS format. The HESSFITS environment has to be set correctly to use this script
  19 + Please see the internal documentation for this purpose:
  20 + https://hess-confluence.desy.de/confluence/display/HESS/HESS+Open+Source+Tools+-+HESS+data+and+IRFs+in+FITS+format
  21 + """
  22 + def __init__(self, *argv):
  23 + """
  24 + Constructor.
  25 + """
  26 + # Set name
  27 + self.name = "cshessobs"
  28 + self.version = "0.1.0"
  29 +
  30 + # Initialise some members
  31 + self.datapath = ""
  32 + self.inmodels = None
  33 + try:
  34 + self.datapath = os.environ["HESSFITS"]
  35 + except KeyError:
  36 + raise gammalib.GException.env_not_found(self.name, "HESSFITS",
  37 + "Please set the environment variable HESSFITS to use this script")
  38 +
  39 + self.outobs = "obs.xml"
  40 + self.nodes = 1
  41 +
  42 + # Make sure that parfile exists
  43 + file = self.parfile()
  44 +
  45 + # Initialise application
  46 + if len(argv) == 0:
  47 + ctools.cscript.__init__(self, self.name, self.version)
  48 + elif len(argv) ==1:
  49 + ctools.cscript.__init__(self, self.name, self.version, *argv)
  50 + else:
  51 + raise TypeError("Invalid number of arguments given.")
  52 +
  53 + # Set logger properties
  54 + self.log_header()
  55 + self.log.date(True)
  56 +
  57 + # Return
  58 + return
  59 +
  60 + def __del__(self):
  61 + """
  62 + Destructor.
  63 + """
  64 + # Write separator into logger
  65 + if self.logTerse():
  66 + self.log("\n")
  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 +
  77 + # Set parfile name
  78 + parfile = self.name+".par"
  79 +
  80 + try:
  81 + pars = gammalib.GApplicationPars(parfile)
  82 + except:
  83 + # Signal that parfile was not found
  84 + print("Parfile \""+parfile+"\" not found. Create default parfile.")
  85 +
  86 + # Check for folders and possible/available configurations
  87 + chains = os.listdir(self.datapath)
  88 +
  89 + # Initialise folder dictionary
  90 + folders = {}
  91 +
  92 + # Loop over possible analysis chains
  93 + for chain in chains:
  94 +
  95 + # Omit hidden files/folders
  96 + if chain[0]=="." or not os.path.isdir(self.datapath+"/"+chain):
  97 + continue
  98 +
  99 + # Initialise dictionary per chain
  100 + folders[chain] = {}
  101 +
  102 + # Read available DST versions
  103 + dstversions = os.listdir(self.datapath+"/"+chain)
  104 +
  105 + # Loop over DST versions
  106 + for dstver in dstversions:
  107 +
  108 + # Omit hidden files/folders
  109 + if dstver[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver):
  110 + continue
  111 +
  112 + # Initialise available cut configurations
  113 + folders[chain][dstver] = []
  114 + configs = os.listdir(self.datapath+"/"+chain+"/"+dstver)
  115 +
  116 + # Loop over configurations
  117 + for config in configs:
  118 +
  119 + # Omit hidden files/folders
  120 + if config[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver+"/"+config):
  121 + continue
  122 +
  123 + # Append configuration to corresponding dictionary
  124 + folders[chain][dstver].append(config)
  125 +
  126 + # Check for appropriate folder structure
  127 + checked_folders = self.check_structure(folders)
  128 + if not len(checked_folders):
  129 + raise KeyError("Appropriate folder structure not found in $HESSFITS. The correct structure can be found on HESS internal documentation")
  130 +
  131 + # Create default parfile
  132 + pars = gammalib.GApplicationPars()
  133 +
  134 + # Create parameters for each available config
  135 + avchains = ""
  136 + for chain in checked_folders:
  137 + avchains += chain +"|"
  138 + avdsts = ""
  139 + for dstver in checked_folders[chain]:
  140 + avdsts+=dstver+"|"
  141 + avconfigs = ""
  142 + for config in checked_folders[chain][dstver]:
  143 + avconfigs += config+"|"
  144 + pars.append(gammalib.GApplicationPar(chain+"_"+dstver+"_config","s","a",checked_folders[chain][dstver][0],avconfigs[:-1],"","Analysis configuration"))
  145 + pars.append(gammalib.GApplicationPar(chain+"_dstver","s","a",checked_folders[chain].keys()[0],avdsts[:-1],"","DST version"))
  146 + pars.append(gammalib.GApplicationPar("chain","s","a",checked_folders.keys()[0],avchains[:-1],"","Analysis chain"))
  147 +
  148 + pars.append(gammalib.GApplicationPar("runlistfile","f","a","runlist.lis","","","Runlist file"))
  149 + pars.append(gammalib.GApplicationPar("inmodel","f","h","NONE","","","Input model XML file (optional)"))
  150 + pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile"))
  151 + pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file"))
  152 + pars.append(gammalib.GApplicationPar("datapath","s","h",self.datapath,"","","Data directoy"))
  153 + pars.append(gammalib.GApplicationPar("nodes","i","h","1","","","Number of free parameters per background model"))
  154 + pars.append(gammalib.GApplicationPar("usetrig","b","h","yes","yes|no","","Pre-estimate background normalisation? Recommended for binned analysis"))
  155 + pars.append_standard()
  156 + pars.save(parfile)
  157 +
  158 + # Return
  159 + return
  160 +
  161 + def get_parameters(self):
  162 + """
  163 + Get parameters from parfile and setup the observation.
  164 + """
  165 + # Get parameters
  166 + self.runlistfile = self["runlistfile"].filename()
  167 +
  168 + # Read models if user input
  169 + self.inmodel = self["inmodel"].filename()
  170 + if not (self.inmodel == "NONE" or self.inmodel == ""):
  171 + self.inmodels = gammalib.GModels(self.inmodel)
  172 +
  173 + # Observation outfile
  174 + self.outobs = self["outobs"].filename()
  175 +
  176 + # Output model file
  177 + self.outmodel = self["outmodel"].filename()
  178 +
  179 + # Read hidden parameters
  180 + self.nodes = self["nodes"].integer()
  181 + self.usetrig = self["usetrig"].boolean()
  182 + self.datapath = self["datapath"].string()
  183 +
  184 + # Read run configuration parameters
  185 + self.chain = self["chain"].string()
  186 + self.dstver = self[self.chain+"_dstver"].string()
  187 + self.config = self[self.chain+"_"+self.dstver+"_config"].string()
  188 +
  189 + # Build directory name prefix
  190 + self.directory = self.datapath+"/"
  191 + self.directory += self.chain +"/"
  192 + self.directory += self.dstver +"/"
  193 + self.directory += self.config +"/"
  194 +
  195 + if not os.path.isdir(self.directory):
  196 + raise gammalib.Exception.runtime_error("Folder structure not found. Please check your HESSFITS environment variable")
  197 +
  198 + self.m_log = False # Logging in client tools
  199 + self.m_debug = False # Debugging in client tools
  200 + self.m_clobber = self["clobber"].boolean()
  201 +
  202 + # Return
  203 + return
  204 +
  205 + def check_structure(self,folders):
  206 +
  207 + # Initialise valid folders
  208 + newfolders = {}
  209 +
  210 + # Check if we have found possible chains
  211 + if len(folders):
  212 +
  213 + # Loop over chains
  214 + for chain in folders:
  215 +
  216 + # Loop over dst versions
  217 + for dst in folders[chain]:
  218 +
  219 + # Loop over configs:
  220 + for config in folders[chain][dst]:
  221 +
  222 + # check if config is valid
  223 + path = self.datapath+"/"+chain+"/"+dst+"/"+config
  224 + if os.path.isdir(path) and len(glob.glob(path+"/run*")):
  225 + if not newfolders.has_key(chain):
  226 + newfolders[chain] = {}
  227 + if not newfolders[chain].has_key(dst):
  228 + newfolders[chain][dst] = []
  229 + newfolders[chain][dst].append(config)
  230 +
  231 + # Return valid configs
  232 + return newfolders
  233 +
  234 + # Create instrument background model
  235 + def hess_inst_background(self,run,aefffile,trgrate=0.0):
  236 +
  237 + # Initialise return spectrum
  238 + spec = 0
  239 +
  240 + # Initialise default value
  241 + value = 1.0
  242 +
  243 + if self.usetrig and trgrate > 0.0:
  244 + # Use value from lookup
  245 + value = 0.00451521513379 * trgrate + 0.162316854264
  246 +
  247 + if self.nodes <= 1:
  248 +
  249 + # Use a constant scaling function
  250 + spec = gammalib.GModelSpectralConst()
  251 + spec["Value"].min(0.1)
  252 + spec["Value"].max(10)
  253 + spec["Value"].value(value)
  254 + if self.nodes < 1:
  255 + spec["Value"].fix()
  256 + else:
  257 + spec["Value"].free()
  258 +
  259 + elif self.nodes == 2:
  260 +
  261 + # User a power law scaling function
  262 + e= gammalib.GEnergy(1.0,"TeV")
  263 + spec = gammalib.GModelSpectralPlaw(1.0,0.0,e)
  264 + spec[0].free()
  265 + spec[0].min(0.01)
  266 + spec[0].max(10.0)
  267 + spec[0].scale(1.0)
  268 + spec[0].value(value)
  269 + spec[1].free()
  270 + spec[1].min(-2.0)
  271 + spec[1].max(2.0)
  272 + spec[1].scale(1.0)
  273 + else:
  274 +
  275 + # Use several energy-dependent scaling factors
  276 + fits = gammalib.GFits(aefffile)
  277 + emin = fits["EFFECTIVE AREA"].real("LO_THRES")
  278 + emax = fits["EFFECTIVE AREA"].real("HI_THRES")
  279 + fits.close()
  280 + spec = gammalib.GModelSpectralNodes()
  281 + bounds = gammalib.GEbounds(self.nodes,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True)
  282 + for i in range(bounds.size()):
  283 + spec.append(bounds.elogmean(i),value)
  284 + for par in spec:
  285 + if "Energy" in par.name():
  286 + par.fix()
  287 + elif "Intensity" in par.name():
  288 + par.min(0.01)
  289 + par.max(10.0)
  290 +
  291 + # Create background model instance
  292 + bck = gammalib.GCTAModelIrfBackground(spec)
  293 +
  294 + # Assign specific run id
  295 + bck.ids(str(int(run)))
  296 +
  297 + # Assign instrument
  298 + bck.instruments("HESS")
  299 +
  300 + # Set name (arbitrary)
  301 + bck.name("bkg_"+run)
  302 +
  303 + # Turn off TS calculation for background model
  304 + bck.tscalc(False)
  305 +
  306 + # Return model
  307 + return bck
  308 +
  309 +
  310 + def execute(self):
  311 + """
  312 + Execute the script.
  313 + """
  314 + # Run the script
  315 + self.run()
  316 +
  317 + # Save residual map
  318 + self.save(self.outobs, self.m_clobber)
  319 +
  320 + # Return
  321 + return
  322 +
  323 + # Get runfolder from runnumber (subject to defined folder structure)
  324 + def get_runfolder(self,runnr):
  325 +
  326 + start = runnr[:-2]+"00"
  327 + if not int(runnr[-3]) % 2 == 0:
  328 + start = str(int(start)-100)
  329 + if len(start) == 5:
  330 + start = "0"+start
  331 + int_end = int(start)+199
  332 + end = str(int_end)
  333 + if len(end) == 5:
  334 + end = "0"+ end
  335 +
  336 + runfold = runnr
  337 + if len(runfold) == 5:
  338 + runfold = "0"+runfold
  339 +
  340 + return self.directory+"/run"+start+"-"+end+"/run"+runfold+"/"
  341 +
  342 + def run(self):
  343 + """
  344 + Run the script.
  345 + """
  346 + # Switch screen logging on in debug mode
  347 + if self.logDebug():
  348 + self.log.cout(True)
  349 +
  350 + # Get parameters
  351 + self.get_parameters()
  352 +
  353 + # Write input parameters into logger
  354 + if self.logTerse():
  355 + self.log_parameters()
  356 + self.log("\n")
  357 +
  358 + # Initialise empty models
  359 + self.models = gammalib.GModels()
  360 +
  361 + # Initialise empty xml file and append an observation list
  362 + self.xml = gammalib.GXml()
  363 + self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\""));
  364 + lib = self.xml.element("observation_list", 0);
  365 +
  366 + # Read runlist from file
  367 + runlist = []
  368 + runfile = open(self.runlistfile)
  369 + for line in runfile.readlines():
  370 + if len(line) == 0:
  371 + continue
  372 + if line[0] == '#':
  373 + continue
  374 + if len(line.split()) > 0:
  375 + runlist.append(line.split()[0])
  376 + else:
  377 + runlist.append('')
  378 + runfile.close()
  379 +
  380 + # Loop over runs
  381 + # Log output
  382 + if self.logTerse():
  383 + self.log("\n")
  384 + self.log.header1("Looping over runs")
  385 +
  386 + for runnr in runlist:
  387 +
  388 + # Unify run number string length
  389 + run = runnr
  390 + if len(run) == 5:
  391 + run = "0"+run
  392 +
  393 + # Get folder in which run files should be located
  394 + runfolder = self.get_runfolder(run)
  395 +
  396 + # construct filenames
  397 + eventfile = runfolder + "hess_events_"+run+".fits.gz"
  398 + aefffile = runfolder + "hess_aeff_2d_"+run+".fits.gz"
  399 + psffile = runfolder + "hess_psf_king_"+run+".fits.gz"
  400 + edispfile = ""#runfolder + "hess_edisp_2d_"+run+".fits.gz"
  401 + bgfile = runfolder + "hess_bkg_offruns_"+run+".fits.gz"
  402 +
  403 + # Check for existence of files
  404 + msg = ""
  405 + skip = False
  406 + if not os.path.isdir(runfolder):
  407 + msg = "Run "+str(int(run))+" does not exist in FITS format"
  408 + skip = True
  409 + elif not os.path.isfile(eventfile):
  410 + msg = "Run "+str(int(run))+" has no eventfile - Run is skipped"
  411 + skip = True
  412 + elif not os.path.isfile(aefffile):
  413 + msg = "Run "+str(int(run))+" has no effective area - Run is skipped"
  414 + skip = True
  415 + elif not os.path.isfile(psffile):
  416 + msg = "Run "+str(int(run))+" has no PSF - Run is skipped"
  417 + skip = True
  418 + elif not os.path.isfile(bgfile):
  419 + msg = "Run "+str(int(run))+" has no background - Run is skipped"
  420 + skip = True
  421 + else:
  422 + msg = "Adding run "+str(int(run))
  423 +
  424 + # Log ouput and skip run if not available
  425 + if self.logTerse():
  426 + self.log(msg+"\n")
  427 + if skip:
  428 + continue
  429 + if self.logExplicit():
  430 + self.log(" Event file: "+eventfile+"\n")
  431 + self.log(" Effective area: "+aefffile+"\n")
  432 + self.log(" Point spread function: "+psffile+"\n")
  433 + self.log(" Energy dispersion: "+edispfile+"\n")
  434 + self.log(" Background file: "+bgfile+"\n")
  435 +
  436 + # Open fits file to determine the observation name
  437 + fits = gammalib.GFits(eventfile)
  438 + object_name = fits["EVENTS"].string("OBJECT")
  439 + trgrate = fits["EVENTS"].real("ZTRGRATE")
  440 + fits.close()
  441 +
  442 + # Append observation to XML and set attributes
  443 + obs = lib.append("observation");
  444 + obs.attribute("name", object_name);
  445 + obs.attribute("id", str(int(run)));
  446 + obs.attribute("instrument", "HESS");
  447 +
  448 + # Append event file
  449 + par = gammalib.GXmlElement("parameter name=\"EventList\"")
  450 + par.attribute("file", eventfile);
  451 + obs.append(par)
  452 +
  453 + # Append effective area
  454 + aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"")
  455 + aeff.attribute("file",aefffile)
  456 +
  457 + # Append PSF
  458 + psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"")
  459 + psf.attribute("file",psffile)
  460 +
  461 + # Append energy dispersion
  462 + edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"")
  463 + edisp.attribute("file",edispfile)
  464 +
  465 + # Append background
  466 + bck = gammalib.GXmlElement("parameter name=\"Background\"")
  467 + bck.attribute("file",bgfile)
  468 +
  469 + # assign events and IRFs to observation
  470 + obs.append(aeff)
  471 + obs.append(psf)
  472 + obs.append(edisp)
  473 + obs.append(bck)
  474 +
  475 + # Append instrumental background model
  476 + self.models.append(self.hess_inst_background(run,aefffile,trgrate))
  477 +
  478 + # Append models provided by 'inmodels' of necessary
  479 + if not self.inmodels == None:
  480 + if self.logTerse():
  481 + self.log("\n")
  482 + self.log.header1("Appending models")
  483 +
  484 + # Loop over input models
  485 + for model in self.inmodels:
  486 + if self.logTerse():
  487 + self.log.header3("Adding model \""+model.name()+"\"")
  488 + self.models.append(model)
  489 +
  490 + # Return
  491 + return
  492 +
  493 + def save(self,outfile,clobber):
  494 +
  495 + # Save observation XML file
  496 + self.xml.save(outfile)
  497 +
  498 + # Save model XML file
  499 + self.models.save(self.outmodel)
  500 +
  501 +# ======================== #
  502 +# Main routine entry point #
  503 +# ======================== #
  504 +if __name__ == '__main__':
  505 + """
  506 + Generates observation file.
  507 + """
  508 + # Create instance of application
  509 + app = cshessobs(sys.argv)
  510 +
  511 + # Open logfile
  512 + app.logFileOpen()
  513 +
  514 + # Execute application
  515 + app.execute()
  516 +
... ...