Thierry Reposeur / ctools

Commit a32545ccc9ab69f5e44f134516480bec31c5e4a2

Authored by Michael Mayer
1 parent b9a29e18

changed import style in csresmap and added a test to test_cscripts.sh.in

Inline Side-by-side
Showing 2 changed files with 85 additions and 37 deletions
scripts/csresmap.py
  View file @a32545c
... ... @@ -19,18 +19,17 @@
19 19 # along with this program. If not, see <http://www.gnu.org/licenses/>.
20 20 #
21 21 # ==========================================================================
22   -from ctools import *
23   -from gammalib import *
  22 +import ctools
  23 +import gammalib
24 24 import sys
25   -import csv
26   -import math
27 25 import tempfile
  26 +import os
28 27  
29 28  
30 29 # ============== #
31 30 # cstsdist class #
32 31 # ============== #
33   -class csresmap(GApplication):
  32 +class csresmap(gammalib.GApplication):
34 33 """
35 34 This class implements the creation of a residual map. It derives
36 35 from the GammaLib::GApplication class which provides support for parameter
... ... @@ -50,19 +49,21 @@ class csresmap(GApplication):
50 49 self.algorithm = "SUB"
51 50  
52 51 # Initialise some members
53   - if isinstance(argv[0],GObservations):
  52 + if isinstance(argv[0],gammalib.GObservations):
54 53 self.obs = argv[0]
  54 + argv = argv[1:]
55 55 else:
56   - self.obs = GObservations()
  56 + self.obs = gammalib.GObservations()
  57 + self.obs.clear()
57 58 self.outfile = ""
58 59 # Make sure that parfile exists
59 60 file = self.parfile()
60 61  
61 62 # Initialise application
62 63 if len(argv) == 0:
63   - GApplication.__init__(self, self.name, self.version)
  64 + gammalib.GApplication.__init__(self, self.name, self.version)
64 65 elif len(argv) ==1:
65   - GApplication.__init__(self, self.name, self.version, argv[1:])
  66 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
66 67 else:
67 68 raise TypeError("Invalid number of arguments given.")
68 69  
... ... @@ -93,32 +94,32 @@ class csresmap(GApplication):
93 94 parfile = self.name+".par"
94 95  
95 96 try:
96   - pars = GApplicationPars(parfile)
  97 + pars = gammalib.GApplicationPars(parfile)
97 98 except:
98 99 # Signal if parfile was not found
99 100 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
100 101  
101 102 # 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"))
  103 + pars = gammalib.GApplicationPars()
  104 + pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition"))
  105 + pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name"))
  106 + pars.append(gammalib.GApplicationPar("outfile","f","a","resmap.fits","","","Output file name"))
  107 + pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
  108 + pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
  109 + pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  110 + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  111 + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)"))
  112 + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)"))
  113 + pars.append(gammalib.GApplicationPar("enumbins","i","h","20","","","Number of energy bins"))
  114 + pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
  115 + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
  116 + pars.append(gammalib.GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"))
  117 + pars.append(gammalib.GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)"))
  118 + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)"))
  119 + pars.append(gammalib.GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels"))
  120 + pars.append(gammalib.GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels"))
  121 + pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)"))
  122 + pars.append(gammalib.GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm"))
122 123 pars.append_standard()
123 124 pars.save(parfile)
124 125  
... ... @@ -133,9 +134,27 @@ class csresmap(GApplication):
133 134  
134 135 # Set observation if not done before
135 136 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()
  137 +
  138 + obsfile = self["infile"].filename()
  139 + try:
  140 +
  141 + self.obs = gammalib.GObservations(obsfile)
  142 + except:
  143 + self.obs.clear()
  144 + observation = gammalib.GCTAObservation()
  145 + observation.load(obsfile)
  146 +
  147 + self.m_irf = self["irf"].string()
  148 + self.m_caldb = self["caldb"].string()
  149 +
  150 + caldb = gammalib.GCaldb()
  151 + if os.path.isdir(self.m_caldb):
  152 + caldb.rootdir(self.m_caldb)
  153 + else:
  154 + caldb.open("cta",self. m_caldb)
  155 +
  156 + observation.response(self.m_irf, caldb);
  157 + self.obs.append(observation)
139 158  
140 159 # Check for models in the container
141 160 # read models from file if there were none
... ... @@ -218,7 +237,7 @@ class csresmap(GApplication):
218 237 self.log.header1("Generate binned map")
219 238  
220 239 # Create countsmap
221   - bin = ctbin(self.obs)
  240 + bin = ctools.ctbin(self.obs)
222 241 bin["nxpix"].integer(self.m_nxpix)
223 242 bin["nypix"].integer(self.m_nypix)
224 243 bin["proj"].string(self.m_proj)
... ... @@ -244,7 +263,7 @@ class csresmap(GApplication):
244 263 self.log.header1("Generate model map")
245 264  
246 265 # create model map
247   - model = ctmodel(self.obs)
  266 + model = ctools.ctmodel(self.obs)
248 267 model["infile"].filename(bin["outfile"].filename())
249 268 model["clobber"].boolean(self.m_clobber)
250 269 if not self.m_modfile == "NONE":
... ... @@ -255,11 +274,11 @@ class csresmap(GApplication):
255 274 model.execute()
256 275  
257 276 # Get model map into GSkymap object
258   - modelmap = GSkymap(model["outfile"].filename())
  277 + modelmap = gammalib.GSkymap(model["outfile"].filename())
259 278 modelmap.stack_maps()
260 279  
261 280 # Initialise resiudalmap from countmap
262   - residualmap = GSkymap(bin["outfile"].filename())
  281 + residualmap = gammalib.GSkymap(bin["outfile"].filename())
263 282 residualmap.stack_maps()
264 283  
265 284  
... ...
test/test_cscripts.sh.in
  View file @a32545c
... ... @@ -131,3 +131,32 @@ else
131 131 exit 1
132 132 fi
133 133 $ECHO " ok"
  134 +
  135 +#
  136 +# Test csresmap
  137 +# =============
  138 +$ECHO @ECHO_N@ "Test csresmap: "
  139 +chmod +x @abs_top_srcdir@/scripts/csresmap.py
  140 +@abs_top_srcdir@/scripts/csresmap.py \
  141 + infile="data/crab_events.fits" \
  142 + outfile="residualmap.fits" \
  143 + srcmdl="data/crab.xml" \
  144 + caldb="irf" \
  145 + irf="cta_dummy_irf" \
  146 + nxpix=50 \
  147 + nypix=50 \
  148 + binsz=0.02 \
  149 + coordsys="CEL" \
  150 + xref=83.63 \
  151 + yref=22.01 \
  152 + proj="CAR" \
  153 + binsz=0.05 > csresmap.log
  154 +$ECHO @ECHO_N@ "."
  155 +if [ -s "residualmap.fits" ]
  156 +then
  157 + $ECHO @ECHO_N@ "."
  158 +else
  159 + $ECHO " residualmap.fits file is not valid"
  160 + exit 1
  161 +fi
  162 +$ECHO " ok"
... ...