#! /usr/bin/env python # ========================================================================== # Residual map generation script. # # Copyright (C) 2014-2016 Michael Mayer # # 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 ctools import sys # ============== # # cstsdist class # # ============== # class csresmap(ctools.cscript): """ This class implements the creation of a residual map. It derives from the ctools.cscript 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 = "csresmap" self.version = "1.1.0" # Initialise some members self.obs = None self.algorithm = "SUB" self.resmap = None self.m_modcube = "None" self.m_outmap = "None" self.m_use_maps = False self.m_skip_binning = False self.m_publish = False # Initialise some members if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): self.obs = argv[0] argv = argv[1:] else: self.obs = gammalib.GObservations() self.obs.clear() self.outfile = "" # Make sure that parfile exists file = self.parfile() # Initialise application if len(argv) == 0: ctools.cscript.__init__(self, self.name, self.version) elif len(argv) ==1: ctools.cscript.__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("inobs","f","a","events.fits","","","Input event list, counts cube, or observation definition XML file")) pars.append(gammalib.GApplicationPar("modcube","f","a","NONE","","","Input model cube file (generated with ctmodel)")) pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Input background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output residual map file")) pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm")) pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)")) pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)")) pars.append(gammalib.GApplicationPar("enumbins","i","h","20","","","Number of energy bins")) pars.append(gammalib.GApplicationPar("ebinfile","f","h","NONE","","","Name of the file containing the energy bin definition")) pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System")) pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method")) pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)")) pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)")) pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels")) pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels")) pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)")) pars.append(gammalib.GApplicationPar("algorithm","s","a","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residual map computation algorithm")) pars.append(gammalib.GApplicationPar("publish","b","h","no","","","Publish residual map on VO Hub?")) pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","csresmap.log","","","Log filename")) pars.save(parfile) # Return return def get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Initialise some flags self.m_use_maps = False self.m_skip_binning = False # First check if the inobs parameter is a counts cube if self.obs.size() == 0 and self["inobs"].filename() != "NONE": filename = gammalib.GFilename(self["inobs"].filename()) if filename.is_fits(): cta = gammalib.GCTAObservation() cta.load(filename) if cta.eventtype() == "CountsCube": self.m_skip_binning = True # If we have a counts cube, then ask whether we also have a model if self.m_skip_binning: self.m_modcube = self["modcube"].filename() if self.m_modcube != "NONE": self.m_use_maps = True # If not two maps are given, proceed to set up observation if not self.m_use_maps: # Set observation if not done before if self.obs.size() == 0: self.require_inobs("csresmap.get_parameters()") self.obs = self.get_observations() # Check if we have exactly one binned CTA observation if self.obs.size() == 1: if self.obs[0].classname() == "GCTAObservation": if self.obs[0].eventtype() == "CountsCube": # Skip ctbin step later on self.m_skip_binning = True # Set models if we have none if self.obs.models().size() == 0: self.obs.models(self["inmodel"].filename()) # Skip query for spatial parameters if a binning is provided in the observation if not self.m_skip_binning: # Read other parameters self.m_xref = self["xref"].real() self.m_yref = self["yref"].real() self.m_emin = self["emin"].real() self.m_emax = self["emax"].real() self.m_enumbins = self["enumbins"].integer() self.m_ebinalg = self["ebinalg"].string() self.m_coordsys = self["coordsys"].string() self.m_proj = self["proj"].string() self.m_nxpix = self["nxpix"].integer() self.m_nypix = self["nypix"].integer() self.m_binsz = self["binsz"].real() # Read energy dispersion flag self.m_edisp = self["edisp"].boolean() # Read algorithm self.m_algorithm = self["algorithm"].string() # Read standard parameters self.m_publish = self["publish"].boolean() self.m_log = False # Logging in client tools self.m_chatter = self["chatter"].integer() self.m_clobber = self["clobber"].boolean() self.m_debug = self["debug"].boolean() # Read ahead output parameters if (self.read_ahead()): self.m_outmap = self["outmap"].filename() # Return return def models(self, models): """ Set model. """ # Copy models self.obs.models(models.clone()) # Return return def execute(self): """ Execute the script. """ # Read ahead output parameters self.read_ahead(True) # Run the script self.run() # Save residual map self.save() # Return return def save(self): """ Save residual map. """ # Write header if self.logTerse(): self.log("\n") self.log.header1("Save residual map") # Continue only filename and residual map are valid if self.m_outmap != "None" and self.resmap != None: # Log file name if self.logTerse(): self.log("Save residual map into \""+self.m_outmap+"\".\n") # Save residual map self.resmap.save(self.m_outmap, self.m_clobber) # Return return def publish(self, name=""): """ Publish residual map. """ # Write header if self.logTerse(): self.log("\n") self.log.header1("Publish sky map") # Continue only if residual map is valid if self.resmap != None: # Set default name is user name is empty if not name: user_name = self.name else: user_name = name # Log map name if self.logTerse(): self.log("Publish residual map \""+user_name+"\".\n") # Publish map self.resmap.publish(user_name) # 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 observation into logger if self.logTerse(): self.log("\n") self.log.header1("Observation") self.log(str(self.obs)) self.log("\n") # Use input file directly if given if self.m_use_maps: countmap = gammalib.GSkyMap(self["inobs"].filename()) modelmap = gammalib.GSkyMap(self.m_modcube) else: if self.m_skip_binning: cta_counts_cube = gammalib.GCTAEventCube(self.obs[0].events().clone()) else: # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate binned map (ctbin)") # Create countsmap bin = ctools.ctbin(self.obs) bin["nxpix"].integer(self.m_nxpix) bin["nypix"].integer(self.m_nypix) bin["proj"].string(self.m_proj) bin["coordsys"].string(self.m_coordsys) bin["xref"].real(self.m_xref) bin["yref"].real(self.m_yref) bin["enumbins"].integer(self.m_enumbins) bin["ebinalg"].string(self.m_ebinalg) bin["emin"].real(self.m_emin) bin["emax"].real(self.m_emax) bin["binsz"].real(self.m_binsz) bin["chatter"].integer(self.m_chatter) bin["clobber"].boolean(self.m_clobber) bin["debug"].boolean(self.m_debug) bin.run() # Retrieve counts cube cta_counts_cube = bin.cube() # Assign GCTAEventCube to skymap countmap = cta_counts_cube.map() # Write header if self.logTerse(): self.log("\n") self.log.header1("Generate model map (ctmodel)") # Create model map model = ctools.ctmodel(self.obs) model.cube(cta_counts_cube) model["chatter"].integer(self.m_chatter) model["clobber"].boolean(self.m_clobber) model["debug"].boolean(self.m_debug) model["edisp"].boolean(self.m_edisp) model.run() # Get model map into GSkyMap object modelmap = model.cube().map().copy() # Store counts map as residual map. Note that we need a # special construct here to avoid memory leaks. This seems # to be a SWIG feature as SWIG creates a new object when # calling bin.cube() #residualmap = bin.cube().map() self.resmap = countmap.copy() self.resmap.stack_maps() modelmap.stack_maps() # Continue calculations depending on given algorithm if self.m_algorithm == "SUB": # Subtract maps self.resmap -= modelmap elif self.m_algorithm == "SUBDIV": # Subtract and divide by model map self.resmap -= modelmap self.resmap /= modelmap #for pixel in modelmap: # if pixel != 0.0: # pixel = 1.0/pixel #self.resmap *= modelmap elif self.m_algorithm == "SUBDIVSQRT": # subtract and divide by sqrt of model map self.resmap -= modelmap self.resmap /= modelmap.sqrt() #for pixel in modelmap: # if pixel != 0.0: # pixel = 1.0/math.sqrt(pixel) #self.resmap *= modelmap else: # Raise error if algorithm is unkown raise TypeError("Algorithm \""+self.m_algorithm+"\" not known") # Optionally publish map if self.m_publish: self.publish() # Return return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': """ Generates residual count map. """ # Create instance of application app = csresmap(sys.argv) # Open logfile app.logFileOpen() # Execute application app.execute()