From 809a67bde8769453016d5c3802b914ad16acb5a1 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 11 Feb 2014 14:42:07 +0100 Subject: [PATCH] Add lsmodel script. --- scripts/Makefile.am | 7 +++++-- scripts/lsmodel.py | 203 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 208 insertions(+), 2 deletions(-) create mode 100755 scripts/lsmodel.py diff --git a/scripts/Makefile.am b/scripts/Makefile.am index 7360197..810837b 100644 --- a/scripts/Makefile.am +++ b/scripts/Makefile.am @@ -4,10 +4,13 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py # Scripts to be distributed in pythondir -pythonscripts = cspull.py cssens.py cstsdist.py obsutils.py +pythonscripts = cspull.py cssens.py cstsdist.py \ + obsutils.py \ + lsmodel.py # Scripts for which symbolic links will be created in bindir -linkscripts = cspull cssens cstsdist +linkscripts = cspull cssens cstsdist \ + lsmodel # Distribute scripts EXTRA_DIST = $(binscripts) $(pythonscripts) diff --git a/scripts/lsmodel.py b/scripts/lsmodel.py new file mode 100755 index 0000000..54666b3 --- /dev/null +++ b/scripts/lsmodel.py @@ -0,0 +1,203 @@ +#! /usr/bin/env python +# ========================================================================== +# This script generates a model map from a counts map and a source model. +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import sys + + +# ============= # +# lsmodel class # +# ============= # +class lsmodel(gammalib.GApplication): + """ + This class implements the model map generation script. It 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 = "lsmodel" + self.version = "0.1.0" + + # 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 script 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("cntmap","f","a","cntmap.fits","","","Counts map")) + pars.append(gammalib.GApplicationPar("expmap","f","a","binned_expmap.fits","","","Exposure map")) + pars.append(gammalib.GApplicationPar("ltcube","f","a","ltcube.fits","","","Livetime cube")) + pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab_model.xml","","","Source model")) + pars.append(gammalib.GApplicationPar("caldb","s","h","","","","Calibration database")) + pars.append(gammalib.GApplicationPar("irf","s","a","P7REP_SOURCE_V15","","","Instrument response function")) + pars.append(gammalib.GApplicationPar("outfile","f","a","model_cube.fits","","","Output file name")) + pars.append_standard() + pars.save(parfile) + + # Return + return + + def get_parameters(self): + """ + Get parameters from parfile and setup the observation. + """ + # Get parameters + self.m_cntmap = self["cntmap"].filename() + self.m_expmap = self["expmap"].filename() + self.m_ltcube = self["ltcube"].filename() + self.m_srcmdl = self["srcmdl"].filename() + self.m_caldb = self["caldb"].string() + self.m_irf = self["irf"].string() + self.m_outfile = self["outfile"].filename() + + # Set some fixed parameters + self.m_log = False # Logging in client tools + self.m_debug = False # Debugging in client tools + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return + + def run(self): + """ + Run the script. + """ + # Switch screen logging on in debug mode + if self.logDebug(): + self.log.cout(True) + + # Get parameters + self.get_parameters() + + # Write input parameters into logger + if self.logTerse(): + self.log_parameters() + self.log("\n") + + # Write header + if self.logTerse(): + self.log("\n") + self.log.header1("Generate model cube") + + # Set LAT observation + obs = gammalib.GLATObservation() + obs.load_binned(self.m_cntmap, self.m_expmap, self.m_ltcube) + obs.response(self.m_irf, self.m_caldb) + + # Set model (use first model for the moment) + models = gammalib.GModels(self.m_srcmdl) + model = models[0] + + # Get deep copy of GLATEventCube from observation + cube = gammalib.GLATEventCube(obs.events()) + + # Fill cube with model values + for i in range(cube.size()): + event = cube[i] + value = model.eval(event, obs) + event.counts(value) + cube[i] = event + + # Log cube + if self.logTerse(): + self.log(str(cube)) + self.log("\n") + + # Write header + if self.logTerse(): + self.log("\n") + self.log.header1("Save model cube") + + # Save model cube + cube.save(self.m_outfile, self.clobber()) + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + """ + Generates a model map from a counts map and a source model. + """ + # Create instance of application + app = lsmodel(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() + \ No newline at end of file -- libgit2 0.21.4