#!/usr/bin/env python import gammalib import ctools import sys import os # ============== # # csxmlbuilder class # # ============== # class csxmlbuilder(gammalib.GApplication): """ This class derives from the GammaLib::GApplication class which provides support for parameter files, command line arguments, and logging. In that way the Python script behaves just as a regular ctool. """ def __init__(self, *argv): """ Constructor. """ # Set name self.name = "csxmlbuilder" self.version = "0.1.0" self.verbose=True # Initialise some members self.datapath = os.environ["HDATA_PATH"] self.outfile = "obs.xml" self.alt_bounds = [0,20,23,27,30,33,37,40,44,49,53,58,64,72,90] self.az_bounds = [-180,-90,90,180] self.nodes = 1 # Make sure that parfile exists file = self.parfile() # Initialise application if len(argv) == 0: gammalib.GApplication.__init__(self, self.name, self.version) elif len(argv) ==1: gammalib.GApplication.__init__(self, self.name, self.version, *argv) else: raise TypeError("Invalid number of arguments given.") # Set logger properties self.log_header() self.log.date(True) # Return return def __del__(self): """ Destructor. """ # Write separator into logger if self.logTerse(): self.log("\n") # Return return def parfile(self): """ Check if parfile exists. If parfile does not exist then create a default parfile. This kluge avoids shipping the cscript with a parfile. """ # Set parfile name parfile = self.name+".par" try: pars = gammalib.GApplicationPars(parfile) except: # Signal if parfile was not found sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("runlistfile","f","a","runlist.list","","","Runlist file")) pars.append(gammalib.GApplicationPar("outfile","f","a","crab.xml","","","observation outfile")) pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","output model file")) pars.append(gammalib.GApplicationPar("datapath","s","h",os.environ["HDATA_PATH"],"","","Data directoy")) pars.append(gammalib.GApplicationPar("nodes","i","h","1","","","Number of parameters per background model")) pars.append_standard() pars.save(parfile) # Return return def get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get parameters self.runlistfile = self["runlistfile"].filename() self.outfile = self["outfile"].filename() self.datapath = self["datapath"].string() self.nodes = self["nodes"].integer() self.outmodel = self["outmodel"].filename() self.m_log = False # Logging in client tools self.m_debug = False # Debugging in client tools self.m_clobber = self["clobber"].boolean() # Return return def hess_inst_background(self,run,emin=0.1,emax=80): spec = 0 if self.nodes <= 1: spec = gammalib.GModelSpectralConst() spec["Value"].min(0.1) spec["Value"].max(10) if self.nodes <1: spec["Value"].fix() else: spec["Value"].free() elif self.nodes == 2: e= GEnergy(1.0,"TeV") spec = gammalib.GModelSpectralPlaw(1.0,0.0,e) spec[0].free() spec[0].min(0.01) spec[0].max(10.0) spec[0].scale(1.0) spec[1].free() spec[1].min(-2.0) spec[1].max(2.0) spec[1].scale(1.0) else: spec = gammalib.GModelSpectralNodes() bounds = gammalib.GEbounds(self.nodes,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"),True) for i in range(bounds.size()): spec.append(bounds.elogmean(i),1.0) for par in spec: if "Energy" in par.name(): par.fix() elif "Intensity" in par.name(): par.min(0.01) par.max(10.0) bck = gammalib.GCTAModelIrfBackground(spec) bck.ids(run) bck.instruments("HESS") bck.name("bck_"+run) bck.tscalc(False) return bck def warning(self, msg): print "*** WARNING *** "+msg def execute(self): """ Execute the script. """ # Run the script self.run() # Save oufile self.save(self.outfile, self.m_clobber) # Return return def _findbins(self,alt,az): altbin=-1 azbin = -1 for i in range(len(self.alt_bounds)-1): if alt >= self.alt_bounds[i] and alt <= self.alt_bounds[i+1]: altbin=i break if az < -90 or az > 90: azbin = 1 else: azbin = 0 return altbin,azbin def run(self): """ Run the script. """ # Switch screen logging on in debug mode if self.logDebug(): self.log.cout(True) self.models = gammalib.GModels() # Get parameters self.get_parameters() self.xml = gammalib.GXml() self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")); lib = self.xml.element("observation_list", 0); k=0 runlist = gammalib.GCsv(self.runlistfile) for i in range(runlist.nrows()): run = runlist.string(i,0) eventfile = self.datapath+"/events/events_"+run+".fits" aefffile = self.datapath+"/Aeff/Aeff_"+run+".fits" psffile = self.datapath+"/psf/psf_"+run+".fits" edispfile = ""#self.datapath+"/events/events_"+run+".fits" bgfile = "" if not os.path.isfile(eventfile): self.warning("Run "+run+" does not exist") continue fits = gammalib.GFits(eventfile) alt = fits[1].real("ALT_PNT") az = fits[1].real("AZ_PNT") ntels = fits[1].integer("N_TELS") object = fits[1].string("OBJECT") fits.close() if alt == 0.0 or az == 0.0 or ntels<1 or ntels>5: if self.verbose: self.warning("Run "+run+" not valid (Alt="+str(alt)+", Az="+str(az)+", Ntels="+str(ntels)+")") continue altbin,azbin = self._findbins(alt, az) bgfile = self.datapath+"/background/hist_alt"+str(altbin)+"_az"+str(azbin)+".fits" if not os.path.exists(bgfile): if self.verbose: self.warning("Run "+run+" has no corresponding off run (run will be skipped, Altbin="+str(altbin)+", Azbin="+str(azbin)+")") continue fits_aeff = gammalib.GFits(aefffile) fits_bg = gammalib.GFits(bgfile) runemax = fits_aeff[1].real("HI_THRES") runemin = fits_bg[1].real("E_THRES") fits_aeff.close() fits_bg.close() obs = lib.append("observation"); obs.attribute("name", object); obs.attribute("id", run); obs.attribute("instrument", "HESS"); obs.attribute("emin",str(runemin)); obs.attribute("emax",str(runemax)); par = gammalib.GXmlElement("parameter name=\"EventList\"") par.attribute("file", eventfile); obs.append(par) aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"") aeff.attribute("file",aefffile) psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"") psf.attribute("file",psffile) edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"") edisp.attribute("file",edispfile) bck = gammalib.GXmlElement("parameter name=\"Background\"") bck.attribute("file",bgfile) obs.append(aeff) obs.append(psf) obs.append(edisp) obs.append(bck) self.models.append(self.hess_inst_background(run, runemin, runemax)) # Return return def save(self,outfile,clobber): self.xml.save(outfile) self.models.save(self.outmodel) # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': """ Generates xml file """ # Create instance of application app = csxmlbuilder(sys.argv) # Open logfile app.logFileOpen() # Execute application app.execute()