Domenico Tiziani / ctools

Commit b9a29e187fd7a022ea20577790484187b57d2fc7

Authored by Michael Mayer
1 parent 9482abe1

added first version of csresmap

Inline Side-by-side
Showing 2 changed files with 323 additions and 2 deletions
scripts/Makefile.am
  View file @b9a29e1
... ... @@ -5,12 +5,12 @@ 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 \
  8 +pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
9 9 lsmodel.py \
10 10 obsutils.py
11 11  
12 12 # Scripts for which symbolic links will be created in bindir
13   -linkscripts = cspull cssens cstsdist cscaldb \
  13 +linkscripts = cspull cssens cstsdist cscaldb csresmap\
14 14 lsmodel
15 15  
16 16 # Distribute scripts
... ...
scripts/csresmap.py 0 → 100755
  View file @b9a29e1
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates the TS distribution for a particular model based
  4 +# on Monte-Carlo simulations.
  5 +#
  6 +# Copyright (C) 2011-2014 Juergen Knoedlseder
  7 +#
  8 +# This program is free software: you can redistribute it and/or modify
  9 +# it under the terms of the GNU General Public License as published by
  10 +# the Free Software Foundation, either version 3 of the License, or
  11 +# (at your option) any later version.
  12 +#
  13 +# This program is distributed in the hope that it will be useful,
  14 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +# GNU General Public License for more details.
  17 +#
  18 +# You should have received a copy of the GNU General Public License
  19 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +#
  21 +# ==========================================================================
  22 +from ctools import *
  23 +from gammalib import *
  24 +import sys
  25 +import csv
  26 +import math
  27 +import tempfile
  28 +
  29 +
  30 +# ============== #
  31 +# cstsdist class #
  32 +# ============== #
  33 +class csresmap(GApplication):
  34 + """
  35 + This class implements the creation of a residual map. It derives
  36 + from the GammaLib::GApplication class which provides support for parameter
  37 + files, command line arguments, and logging. In that way the Python
  38 + script behaves just as a regular ctool.
  39 + """
  40 + def __init__(self, *argv):
  41 + """
  42 + Constructor.
  43 + """
  44 + # Set name
  45 + self.name = "csresmap"
  46 + self.version = "0.1.0"
  47 +
  48 + # Initialise some members
  49 + self.obs = None
  50 + self.algorithm = "SUB"
  51 +
  52 + # Initialise some members
  53 + if isinstance(argv[0],GObservations):
  54 + self.obs = argv[0]
  55 + else:
  56 + self.obs = GObservations()
  57 + self.outfile = ""
  58 + # Make sure that parfile exists
  59 + file = self.parfile()
  60 +
  61 + # Initialise application
  62 + if len(argv) == 0:
  63 + GApplication.__init__(self, self.name, self.version)
  64 + elif len(argv) ==1:
  65 + GApplication.__init__(self, self.name, self.version, argv[1:])
  66 + else:
  67 + raise TypeError("Invalid number of arguments given.")
  68 +
  69 + # Set logger properties
  70 + self.log_header()
  71 + self.log.date(True)
  72 +
  73 + # Return
  74 + return
  75 +
  76 + def __del__(self):
  77 + """
  78 + Destructor.
  79 + """
  80 + # Write separator into logger
  81 + if self.logTerse():
  82 + self.log("\n")
  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 + # Set parfile name
  93 + parfile = self.name+".par"
  94 +
  95 + try:
  96 + pars = GApplicationPars(parfile)
  97 + except:
  98 + # Signal if parfile was not found
  99 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  100 +
  101 + # Create default parfile
  102 + pars = GApplicationPars()
  103 + pars.append(GApplicationPar("infile","f","h","obs.xml","","","Observation definition"))
  104 + pars.append(GApplicationPar("srcmdl","f","h","crab.xml","","","model file name"))
  105 + pars.append(GApplicationPar("outfile","f","h","resmap.fits","","","Output file name"))
  106 + pars.append(GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
  107 + pars.append(GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
  108 + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  109 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  110 + pars.append(GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)"))
  111 + pars.append(GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)"))
  112 + pars.append(GApplicationPar("enumbins","i","h","20","","","Number of energy bins"))
  113 + pars.append(GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
  114 + pars.append(GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
  115 + pars.append(GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"))
  116 + pars.append(GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)"))
  117 + pars.append(GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)"))
  118 + pars.append(GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels"))
  119 + pars.append(GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels"))
  120 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)"))
  121 + pars.append(GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm"))
  122 + pars.append_standard()
  123 + pars.save(parfile)
  124 +
  125 + # Return
  126 + return
  127 +
  128 + def get_parameters(self):
  129 + """
  130 + Get parameters from parfile and setup the observation.
  131 + """
  132 + # Get parameters
  133 +
  134 + # Set observation if not done before
  135 + if self.obs.size() == 0:
  136 + self.obs = GObservations(self["infile"].filename)
  137 + self.m_irf = self["irf"].string()
  138 + self.m_caldb = self["caldb"].string()
  139 +
  140 + # Check for models in the container
  141 + # read models from file if there were none
  142 + # Note that the models should be optimised for
  143 + # meaningful residual map
  144 + if self.obs.models().size() == 0:
  145 + self.obs.models(self["srcmdl"].filename())
  146 + self.m_outfile = self["outfile"].filename()
  147 +
  148 + self.m_modfile = self["modfile"].filename()
  149 + self.m_binfile = self["binfile"].filename()
  150 +
  151 + self.m_xref = self["xref"].real()
  152 + self.m_yref = self["yref"].real()
  153 + self.m_emin = self["emin"].real()
  154 + self.m_emax = self["emax"].real()
  155 + self.m_enumbins = self["enumbins"].integer()
  156 + self.m_ebinalg = self["ebinalg"].string()
  157 + self.m_coordsys = self["coordsys"].string()
  158 + self.m_proj = self["proj"].string()
  159 + self.m_nxpix = self["nxpix"].integer()
  160 + self.m_nypix = self["nypix"].integer()
  161 + self.m_binsz = self["binsz"].real()
  162 + self.m_algorithm = self["algorithm"].string()
  163 +
  164 + # Set some fixed parameters
  165 + self.m_log = False # Logging in client tools
  166 + self.m_debug = False # Debugging in client tools
  167 + self.m_clobber = self["clobber"].boolean()
  168 +
  169 + # Return
  170 + return
  171 +
  172 + def models(self, models):
  173 + """
  174 + Set model.
  175 + """
  176 + # Copy models
  177 + self.obs.models(models.clone())
  178 +
  179 + # Return
  180 + return
  181 +
  182 + def execute(self):
  183 + """
  184 + Execute the script.
  185 + """
  186 + # Run the script
  187 + self.run()
  188 +
  189 + # Return
  190 + return
  191 +
  192 + def run(self):
  193 + """
  194 + Run the script.
  195 + """
  196 + # Switch screen logging on in debug mode
  197 + if self.logDebug():
  198 + self.log.cout(True)
  199 +
  200 + # Get parameters
  201 + self.get_parameters()
  202 +
  203 + # Write input parameters into logger
  204 + if self.logTerse():
  205 + self.log_parameters()
  206 + self.log("\n")
  207 +
  208 + # Write observation into logger
  209 + if self.logTerse():
  210 + self.log("\n")
  211 + self.log.header1("Observation")
  212 + self.log(str(self.obs))
  213 + self.log("\n")
  214 +
  215 + # Write header
  216 + if self.logTerse():
  217 + self.log("\n")
  218 + self.log.header1("Generate binned map")
  219 +
  220 + # Create countsmap
  221 + bin = ctbin(self.obs)
  222 + bin["nxpix"].integer(self.m_nxpix)
  223 + bin["nypix"].integer(self.m_nypix)
  224 + bin["proj"].string(self.m_proj)
  225 + bin["coordsys"].string(self.m_coordsys)
  226 + bin["xref"].real(self.m_xref)
  227 + bin["yref"].real(self.m_yref)
  228 + bin["enumbins"].integer(self.m_enumbins)
  229 + bin["ebinalg"].string(self.m_ebinalg)
  230 + bin["emin"].real(self.m_emin)
  231 + bin["emax"].real(self.m_emax)
  232 + bin["binsz"].real(self.m_binsz)
  233 + bin["clobber"].boolean(self.m_clobber)
  234 + if not self.m_binfile == "NONE":
  235 + bin["outfile"].filename(self.m_binfile)
  236 + else:
  237 + os,filename = tempfile.mkstemp()
  238 + bin["outfile"].filename(filename)
  239 + bin.execute()
  240 +
  241 + # Write header
  242 + if self.logTerse():
  243 + self.log("\n")
  244 + self.log.header1("Generate model map")
  245 +
  246 + # create model map
  247 + model = ctmodel(self.obs)
  248 + model["infile"].filename(bin["outfile"].filename())
  249 + model["clobber"].boolean(self.m_clobber)
  250 + if not self.m_modfile == "NONE":
  251 + model["outfile"].filename(self.m_modfile)
  252 + else:
  253 + os,filename = tempfile.mkstemp()
  254 + model["outfile"].filename(filename)
  255 + model.execute()
  256 +
  257 + # Get model map into GSkymap object
  258 + modelmap = GSkymap(model["outfile"].filename())
  259 + modelmap.stack_maps()
  260 +
  261 + # Initialise resiudalmap from countmap
  262 + residualmap = GSkymap(bin["outfile"].filename())
  263 + residualmap.stack_maps()
  264 +
  265 +
  266 +
  267 +
  268 + # continue calculations depending on given algorithm
  269 + if self.m_algorithm == "SUB":
  270 +
  271 + # Subtract maps
  272 + residualmap -= modelmap
  273 +
  274 + elif self.m_algorithm == "SUBDIV":
  275 +
  276 + # Subtract and divide by model map
  277 + residualmap -= modelmap
  278 + residualmap /= modelmap
  279 +
  280 + elif self.m_algorithm == "SUBDIVSQRT":
  281 +
  282 + # subtract and divide by sqrt of model map
  283 + residualmap -= modelmap
  284 + residualmap /= modelmap.sqrt()
  285 +
  286 + else:
  287 +
  288 + # Rais error if algorithm is unkown
  289 + raise TypeError("Algorithm \""+self.m_algorithm+"\" not known")
  290 +
  291 + # Save files
  292 + residualmap.save(self.m_outfile,self.m_clobber)
  293 +
  294 + # Remove binned and model files in case they shouldnt be saved
  295 + if self.m_binfile == "NONE":
  296 + os.remove(bin["outfile"].filename())
  297 + if self.m_modfile == "NONE":
  298 + os.remove(model["outfile"].filename())
  299 +
  300 + # Return
  301 + return
  302 +
  303 +
  304 +
  305 +
  306 +# ======================== #
  307 +# Main routine entry point #
  308 +# ======================== #
  309 +if __name__ == '__main__':
  310 + """
  311 + Generates residual count map.
  312 + """
  313 + # Create instance of application
  314 + app = csresmap(sys.argv)
  315 +
  316 + # Open logfile
  317 + app.logFileOpen()
  318 +
  319 + # Execute application
  320 + app.execute()
  321 +
0 322 \ No newline at end of file
... ...